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In  multi-target  tracking  problems  such  as  those  found  in  high-energy 
particle  physics,  fluid  mechanics,  and  ballistic  missile  defense,  the  common  ob¬ 
jective  is  to  separate  the  data  into  observations  associated  with  individual  targets 
and  to  use  this  data  to  estimate  the  targets’  trajectories.  In  defense  related  ap¬ 
plications,  it  is  necessary  to  have  algorithms  which  are  computationally  efficient, 
robust,  and  minimize  data  storage  requirements.  Recently  developed  approaches 
in  the  field  of  multi-target  tracking,  however,  have  been  shown  to  have  significant 
computational  disadvantages. 

In  this  study,  non-hierarchical  clustering  methods  are  combined  with 
computationally  efficient  algorithms  such  as  those  used  to  solve  assignment  and 
quadratic  programming  problems  to  provide  an  integrated  procedure  which  is 
computationally  efficient,  minimizes  data  storage  requirements,  and  gives  a  rea¬ 
sonable  estimate  of  the  number  of  targets.  Combined  with  a  sequential  estimation 
filter  such  as  the  extended  Kalman  filter,  the  procedure  can  provide  estimates  of 
a  target’s  state  and  state  covariance  after  three  observations  and  continuously 
maintain  updated  target  state  estimates  in  real  time. 
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Empirical  results  based  on  100  targets  in  ballistic  trajectories  have 
demonstrated  this  method’s  effectiveness  by  properly  clustering  data  with  four 
measurement  attributes  (range,  range  rate,  azimuth,  and  elevation)  in  over  98 
percent  of  the  cases.  Its  robustness  is  manifested  by  the  fact  that  these  results 
apply  to  scenarios  with  20  percent  missing  data  and  biases  of  up  to  one  arc  minute 
in  the  sensor  attitude  and  0.5  seconds  in  the  sensor  clock.  And  its  capability  to 
track  in  real  time  is  demonstrated  with  a  duty  cycle  of  less  than  five  percent. 
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Chapter  1 


Background 


1.1  Introduction 

Today,  many  applications  require  the  tracking  of  a  large  set  of  uniden¬ 
tified  targets.  Among  these  are  applications  in  high-energy  particle  physics,  fluid 
mechanics,  and  ballistic  missile  defense.  Crucial  to  the  concept  of  a  space-based 
ballistic  missile  defense  under  the  proposed  Strategic  Defense  Initiative  (SDI)  is 
the  development  of  a  system  with  the  capability  to  detect,  classify,  and  predict 
the  motion  of  a  large  number  of  unidentified  targets.  Composed  of  orbiting  sen¬ 
sor  platforms  and  linked  through  the  Command,  Control,  Communications,  and 
Intelligence  (C3I)  element,  this  system  must  not  only  be  able  to  handle  multi¬ 
ple  targets,  but  also  to  integrate  the  combined  data  from  multiple  sensors.  This 
data  must  be  combined  in  such  a  way  as  to  present  a  realistic  picture  of  the 
scenario  underway  so  that  limited  resources  may  be  directed  in  an  appropriate 
response.  Because  the  time  for  such  a  response  in  an  Intercontinental  Ballistic 
Missile  (ICBM)  attack  scenario  is  so  short,  efficient,  robust,  and  accurate  algo¬ 
rithms  capable  of  handling  multi-target,  multi-sensor  data  in  real  time  are  critical 
to  achieving  a  successful  ICBM  defense. 

1.2  Problem  Statement 

Nominally,  this  detection-estimation  system  will  be  comprised  of  two 
elements.  The  first  element  is  a  constellation  of  observation  satellites  placed  in 
orbital  configurations  which  allow  suitable  coverage  of  the  areas  of  interest  (the 
ICBMs’  trajectories  from  their  launch  sites  to  the  anticipated  impact  points).  The 
satellites’  onboard  sensors  must  be  capable  of  providing  time-tagged  observations 
of  a  number  of  attributes  of  the  targets  in  its  field-of-view.  Typically,  these 


attributes  might  include  range  and  range-rate  (for  active  sensors)  and  azimuth 
and  elevation  (for  both  active  and  passive  sensors). 

In  a  typical  ballistic  missile  defense  scenario,  these  orbiting  satellite  sen¬ 
sor  platforms  will  survey  the  earth  from  altitudes  of  several  thousand  kilometers. 
Upon  launch  of  an  ICBM  attack,  each  sensor  may  detect  as  many  as  100  tar¬ 
gets  within  its  field-of-view,  and  these  targets  will  likely  be  closely  spaced  and 
have  similar  attributes.  Not  only  will  there  be  uncertainty  associated  with  the 
measurements  of  the  targets’  attributes  due  to  sensor  limitations,  but  some  ob¬ 
servations  will  be  lost  due  to  spurious  measurements  or  unobservable  conditions 
relating  to  the  sensor-target  geometry.  In  addition,  there  will  be  uncertainties 
associated  with  both  the  sensor  attitude  and  sensor  clock  (position). 

The  second  element  of  the  detection-estimation  system  is  the  C3I  site 
where  the  data  from  various  sensors  is  combined.  This  site  might  be  a  land-based 
command  center  or  a  space-based  battle  station  (perhaps  even  co-located  with 
one  of  the  observation  sensors). 

Between  these  two  system  elements,  four  basic  tasks  must  be  accom¬ 
plished: 

•  Separate  the  available  data  into  tracks  associated  with  individual  targets, 

•  Correlate/ combine  tracks  from  various  sensors, 

•  Estimate  the  targets’  state  at  some  reference  epoch,  and 

•  Predict  (track)  the  targets’  state  at  some  future  epoch. 

The  order  in  which  these  tasks  are  listed  should  not  be  taken  to  imply 
a  sequential  relationship.  In  fact,  how  these  tasks  are  performed  will  determine 
the  overall  complexity  of  the  multi-target,  multi-sensor  tracking  problem.  In 
addition,  just  what  processing  is  done  and  by  which  element  is  a  question  of 
distributed  estimation. 

Further  complicating  an  already  difficult  problem  are  uncertainties  in 
the  numbers  of  targets  visible  at  each  sensor  or  jointly  visible  at  any  subset  of 
sensors,  noisy  and  spurious  measurements,  and  the  nonlinear  measurements  and 
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nonlinear  dynamics  of  ballistic  flight.  Additional  constraints  may  arise  due  to 
the  need  to  decentralize  data  processing  to  enhance  survivability. 

1.3  Multi-Target  Tracking 

As  a  result  of  the  many  complexities  involved  in  the  multi-target  track¬ 
ing  problem,  a  wide  range  of  methods  have  been  developed  in  an  attempt  to 
handle  these  four  tasks  and  their  associated  difficulties.  Research  in  multi-target 
tracking  has  concentrated  in  four  primary  areas: 

•  Track  initiation, 

•  Track  maintenance, 

•  Sensor-to-sensor  correlation,  and 

•  Improved  estimation  methods. 

These  four  areas  correspond  roughly  to  the  four  tasks  performed  by  the  detection- 
estimation  system. 

Methods  for  handling  track  initiation  and  track  maintenance  often  have 
much  in  common  and,  as  such,  have  actually  been  handled  as  different  mani¬ 
festations  of  the  same  problem.  Historically,  approaches  to  this  problem  can  be 
classed  as  Bayesian  or  non-Bayesian. 

Initially,  research  focused  on  what  are  now  known  as  non-Bayesian 
methods.  Led  by  the  pioneering  work  of  Sittler  in  1964,  these  methods  include 
(1)  tracking  via  data  association,  (2)  track-split  filtering,  and  (3)  the  maximum 
likelihood  method.  In  tracking  via  data  association,  Sittler  [47]  devised  a  method 
whereby  whenever  more  than  one  sensor  measurement  was  observed  in  the  neigh¬ 
borhood  of  a  predicted  measurement,  the  current  track  was  split.  Trajectories 
whose  maximum  likelihood  function  fell  below  a  certain  threshold  were  dropped 
from  further  consideration.  This  method  handles  both  track  initiation  and  track 
termination,  as  well  as  false  alarms  and  missing  measurements. 

In  1975,  Smith  and  Buechler  [48]  expanded  on  Sittler’s  approach  within 
the  framework  of  Kalman  filtering  (which  was  not  in  common  use  in  1964)  for  an 
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application  to  radar  tracking.  In  the  track-splitting  filter,  Sittler’s  concept  of  a 
neighborhood  was  now  a  validation  region  which  was  derived  from  the  innovation 
covariance  matrix  obtained  from  the  standard  Kalman  filter.  However,  both 
Sittler’s  method  and  that  of  Smith  and  Buechler  were  considered  impractical 
because  the  exponential  growth  in  the  number  of  trajectories  would  saturate  the 
memory  and  computational  capability  of  even  the  largest  computers. 

Stein  and  Blackman  [52]  further  modernized  Sittler’s  work  in  the  devel¬ 
opment  of  the  maximum  likelihood  method.  They  used  a  suboptimal  sequential 
method  which  selected  only  the  most  likely  assignment  of  targets  and  measure¬ 
ments  from  each  data  set  or  scan ,  thereby  mitigating  the  trajectory  growth  prob¬ 
lem.  Morefield  [39]  extended  this  approach  by  partitioning  the  data  into  mutually 
exclusive  and  exhaustive  sets  of  feasible  tracks  and  formulating  a  0-1  integer  pro¬ 
gram.  While  Morefield’s  work  put  track  initiation  and  maintenance  on  a  more 
solid  theoretical  basis,  its  application  still  has  large  computational  and  memory 
requirements  in  a  dense  target  environment. 

Initial  work  with  Bayesian  approaches  began  with  the  nearest  neighbor 
filter.  Sea  [42],  Singer  and  Stein  [44],  and  Singer  and  Sea  [45]  used  the  nearest 
neighbor  of  a  predicted  measurement  and  modified  the  Kalman  filter  to  account 
for  the  a  priori  probability  that  this  measurement  might  be  spurious  (Bar-Shalom 
[10]).  However,  it  was  discovered  that  this  filter  can  easily  lose  the  target  in  a 
cluttered  environment. 

Work  by  Jaffer  and  Bar-Shalom  [31]  and  Bar-Shalom  and  Jaffer  [6]  led  to 
the  development  of  the  probability  data  association  filter  (PDAF)  by  Bar-Shalom 
[8],  Bar-Shalom  and  Tse  [7,9],  and  Bar-Shalom  and  Birmiwal  [11].  A  suboptimal 
Bayesian  approach,  the  PDAF  sequentially  incorporates  clusters  of  measurements 
into  a  track  by  attaching  to  each  cluster  an  a  posteriori  probability  of  being 
correct.  This  is  important  because  the  standard  formulation  of  the  Kalman  filter 
is  optimal  only  when  there  is  no  possibility  of  incorrect  assignments  being  made 
to  a  track.  As  a  result,  estimates  and  covariances  in  the  PDAF  account  for  the 
measurement  origin  uncertainty  rather  than  being  conditioned  on  the  “accepted” 
tracks  being  true.  The  primary  limitation  of  the  PDAF  is  that  it  only  tracks  a 
single  target  in  a  multiple  target  or  cluttered  environment. 
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More  recently,  the  joint  PDAF  of  Fortmann,  Bar-Shalom,  and  Scheffe 
[27],  Chang  and  Bar-Shalom  [18,19],  and  Chang,  Chong,  and  Bar-Shalom  [20]  is 
used  to  jointly  compute  the  probabilities  for  all  targets  and  measurements  that 
form  a  cluster.  In  it,  posterior  probabilities  of  the  joint  probability  distribution 
function  are  conditioned  on  all  measurements  up  to  the  present,  allowing  multiple 
targets  (resolved  or  unresolved)  to  be  tracked  in  a  cluttered  environment. 

An  optimal  Bayesian  approach  was  developed  by  Singer,  Sea,  and  House- 
wright  [46]  for  a  single  target  in  a  cluttered  environment.  The  major  difference 
between  the  optimal  Bayesian  approach  and  the  PDAF  is  that  decomposition  of 
the  state  estimate  is  accomplished  in  terms  of  all  combinations  of  measurements 
from  initial  to  present  time  rather  than  in  terms  of  the  latest  measurements  only. 
That  is,  the  state  estimate  is  determined  based  upon  all  possible  track  histo¬ 
ries  for  a  given  target.  Again,  the  major  difficulty  with  each  of  these  last  two 
approaches  is  exponential  memory  growth. 

In  a  multiple  sensor  environment,  data  from  the  various  sensors  must 
somehow  be  combined  to  correlate  targets  sets  which  may  be  visible  to  sev¬ 
eral  sensors.  Methods  for  sensor-to-sensor  correlation  attributed  to  Singer  and 
Kanyuck  [43],  Stein  and  Blackman  [52],  Bowman  [14],  and  Chang  and  Youens 
[16],  are  generally  extensions  of  the  maximum  likelihood  approaches  used  for  the 
track  maintenance  problem.  There  are  two  primary  approaches  to  correlating 
target  sets  from  multiple  sensors.  The  first  is  to  map  all  observations  from  all 
sensors  into  a  common  measurement  space  and  apply  any  of  the  tracking  methods 
used  for  the  single-sensor  case.  However,  this  approach  will  work  only  with  mea¬ 
surement  sets  which  permit  unambiguous  mappings  into  a  common  measurement 
space.  The  second  approach  is  to  form  single-sensor  tracks  and  then  correlate  the 
tracks  from  various  sensors  via  pattern  recognition  or  matchings  of  target  state 
estimates. 

Finally,  because  of  the  likelihood  of  imperfect  correlation  of  observations 
with  clusters,  it  is  necessary  to  develop  tracking  filters  which  incorporate  these 
correlation  errors  in  the  update  of  the  error  covariance  matrix.  The  efforts  of 
Singer  and  Stein  [44],  Jaffer  and  Bar-Shalom  [31],  Singer  and  Sea  [45],  and  Singer, 
Sea,  and  Housewright  [46]  resulted  in  the  development  of  filtering  algorithms 
which  use  the  a  posteriori  probability  that  an  observation  originated  from  a 


specific  target  in  a  dense  multi-target  environment.  In  addition,  Tse,  Larson, 
and  Bar-Shalom  [56]  and  Chang  [15]  considered  the  specific  problem  of  estimation 
with  angles-only  measurements.  The  particular  set  of  measurements  available  can 
play  a  significant  role  in  both  choosing  the  most  appropriate  tracking  method 
and  developing  the  tracking  filter  due  to  problems  of  marginal  observability.  The 
angles-only  case  is  particularly  difficult  in  this  respect. 

1.4  Clustering 

In  an  attempt  to  mitigate  some  of  the  computational  complexity  of  the 
algorithms  discussed  in  the  previous  section,  recent  investigators  have  examined 
the  application  of  a  broad  range  of  methods  from  the  classification  field  (Tapley 
et  al.  [53,54],  Balakrishnan  et  al.  [4]).  These  methods  are  collectively  known  as 
clustering  methods. 

The  first  definitive  results  in  clustering  were  produced  by  Sokal  and 
Sneath  [50]  and  refined  by  Sneath  and  Sokal  [49]  in  the  field  of  numerical  taxon¬ 
omy.  While  these  references  are  devoted  primarily  to  the  biological  sciences,  the 
approach  applies  to  all  types  of  clustering.  The  clustering  algorithms  presented 
in  these  and  subsequent  references  by  Anderberg  [2]  and  Romesburg  [41]  all  share 
a  common  framework. 

First,  the  data  to  be  clustered  is  standardized  based  on  some  figure 
of  merit  which  considers  the  relative  importance  of  the  various  types  of  mea¬ 
surements.  Once  this  standardization  is  completed,  a  similarity  or  dissimilarity 
coefficient  is  computed  between  each  pair  of  measurements.  The  resulting  matrix 
is  known  as  a  resemblance  matrix.  There  are  many  methods  for  computing  its 
coefficients.  Typically,  these  coefficients  are  metrics  such  as  Euclidean  or  average 
Euclidean  distance  in  the  measurement  (attribute)  space.  Although  other  types 
of  coefficients  exist,  they  are  not  appropriate  to  this  endeavor  because  of  their 
inability  to  discriminate  the  types  of  clusters  encountered  in  the  multi-target 
tracking  problem  (those  with  additive  or  proportional  translations). 

Given  the  resemblance  matrix,  there  are  two  approaches  to  forming 
clusters:  hierarchical  and  non-hierarchical.  The  earliest  approaches  to  clustering 
involved  hierarchical  clustering,  wherein  clusters  are  “built  up”  from  the  data  un- 


til  one  large  cluster  is  formed.  The  most  common  of  the  many  methods  available 
involve  linking  clusters  through  the  use  of  either  weighted  or  unweighted  pair- 
group  procedures  using  arithmetic  averages  (WPGMA  or  UPGMA).  In  these 
methods,  links  are  developed  based  on  the  smallest  (weighted  or  unweighted) 
average  of  all  the  distances1  between  the  points  in  a  pair  of  clusters. 

Other  common  methods  include  the  single  linkage,  complete  linkage, 
and  centroid  methods.  Single  linkage  builds  clusters  through  association  of  the 
shortest  distance  between  the  closest  points  in  any  two  clusters  while  complete 
linkage  builds  clusters  based  on  association  of  the  shortest  distance  between  the 
farthest  points  in  any  two  clusters.  Single  linkage  may  be  thought  of  as  a  “near¬ 
est  neighbor”  approach  whereas  complete  linkage  is  a  “strongest  association” 
approach.  Centroid  methods  offer  a  compromise  between  these  two  extremes, 
building  linkages  based  on  the  shortest  distance  between  the  centroids  of  existing 
clusters.  Of  the  hierarchical  methods  discussed,  the  single  linkage  method  is  the 
most  appropriate  for  the  multi-target  tracking  problem  because  of  its  tendency 
to  form  clusters  which  are  chains  of  data  points  (a  feature  which  is  normally 
considered  a  drawback  to  this  method). 

The  tree  formed  by  any  of  these  clustering  methods,  which  shows  how 
the  clusters  are  linked  and  at  what  level,  is  known  as  a  phenogram  or  dendro¬ 
gram.  In  Figure  1.1,  the  phenogram  shows  four  objects  and  how  they  are  related. 
The  level  at  which  objects  are  linked  together  represents  their  similarity,  with  the 
lowest  links  indicating  the  strongest  similarities.  To  complete  the  hierarchical  ap¬ 
proach  the  phenogram  must  be  “split”  at  some  level  to  decide  how  many  clusters 
exist,  a  requirement  key  to  the  multi-target  tracking  problem  where  the  number 
of  targets  is  unknown.  Depending  on  what  level  the  phenogram  in  Figure  1.1  is 
split,  will  determine  how  many  distinct  clusters  exist.  Splitting  at  Level  1  yields 
four  clusters  while  splitting  at  Level  2  yields  only  two. 

In  hierarchical  clustering  methods,  a  data  set  of  n  observations  yields  n 
nested  classifications  ranging  from  one  cluster  with  n  observations  to  n  clusters 
with  one  observation  each.  Non-hierarchical  methods,  on  the  other  hand,  cluster 

lrrhe  term  distance  is  used  in  this  discussion  to  imply  a  similarity  or  dissimilarity  coefficient. 
Smaller  distances  refer  to  similar  coefficients  while  larger  distances  refer  to  the  opposite. 
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Figure  1.1:  Sample  Phenogram 

observations  into  k  clusters,  where  k  is  either  specified  a  priori  or  is  determined 
as  part  of  the  method  used  (Anderberg  [2]).  These  methods  enjoy  an  advantage 
over  hierarchical  methods  since  it  is  not  necessary  for  them  to  store  the  similarity 
matrix  or  even  the  data  set  since  the  data  is  typically  processed  serially.  It  is, 
therefore,  possible  to  cluster  much  larger  data  sets  with  non-hierarchical  methods. 

The  majority  of  non-hierarchical  clustering  methods  involve  the  use  of 
seed  points.  There  are  many  ways  of  seeding  clusters.  MacQueen  [36]  suggested 
choosing  the  first  k  observations  as  seeds  while  McRae  [37]  chose  k  random  obser¬ 
vations.  Forgy  [26]  partitioned  the  data  into  k  mutually  exclusive  and  exhaustive 
sets  and  used  the  set  centroids  as  the  seed  points.  A  more  intuitively  appealing 
approach  was  used  by  Astrahan  [3]  wherein  “densities”  were  calculated  for  each 
data  point  and  the  points  with  the  k  highest  “densities”  were  selected  as  seeds. 
In  a  similar  approach,  Ball  and  Hall  [5]  chose  the  first  seed  as  the  centroid  of  the 
data  set.  Additional  seed  points  were  added  while  processing  the  data  if  they 
were  more  than  some  set  distance  from  all  existing  seed  points. 

Once  the  seed  points  have  been  determined,  clusters  are  built  around 
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them.  Forgy’s  method  [26]  and  Jancey’s  variant  [32]  assign  observations  to  the 
closest  seed  while  MacQueen’s  fc-means  method  [36]  assigns  observations  to  the 
cluster  with  the  nearest  centroid.  Usually  the  clustering  is  reapplied  after  gen¬ 
erating  new  seed  points  based  on  the  current  partition,  such  as  the  new  cluster 
centroids,  and  repeated  until  some  convergence  criterion  is  satisfied. 

MacQueen  [36]  and  Wishart  (Anderberg  [2])  developed  methods  which 
permit  variable  numbers  of  clusters.  In  these  methods,  clusters  are  merged  if 
their  seeds  are  within  some  pre-specified  distance  of  each  other.  New  clusters 
are  formed  when  observations  are  found  to  be  beyond  some  (usually  different) 
distance  from  the  existing  cluster  seed  points  or  centroids.  As  with  the  fixed 
number  of  clusters  methods,  the  process  is  repeated  until  convergence. 

Finally,  many  authors  have  developed  methods  which  propose  criteria 
for  evaluating  whether  movements  of  individual  observations  result  in  an  overall 
improvement  of  a  partition  (Anderberg  [2],  Spath  [51]).  These  criteria  are  based 
on  multivariate  statistical  analysis  techniques,  such  as  linear  discriminant  analysis 
and  multivariate  analysis  of  variance.  The  principal  criteria  used  are: 

•  Minimize  trace  W, 

•  Minimize  |W|/|T|  or  maximize  jTJ/|W|, 

•  Maximize  the  largest  eigenvalue  of  W~XB ,  and 

•  Maximize  the  trace  of  W-1  B, 

where  T  is  the  total  scatter  matrix,  W  is  the  within  cluster  scatter  matrix,  and 
B  is  the  between  cluster  scatter  matrix.  It  can  be  shown  that  the  three  matrices 
satisfy  the  relation  T  =  B  +  W  (Anderberg  [2]).  These  criteria  are  generally 
applied  to  the  non-hierarchical  methods  discussed  above  as  tests  of  convergence. 

None  of  the  currently  employed  clustering  methods  are  designed  to  ex¬ 
plicitly  handle  temporal  data.  Historically,  clustering  methods  were  developed  to 
segregate  data  into  distinct  classifications.  Because  the  desire  is  to  group  those 
observations  which  are  most  similar,  the  algorithms  in  use  tend  to  generate  hyper- 
spherical  clusters  in  the  attribute  space.  As  applied  to  the  multi-target  tracking 
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problem,  however,  this  tendency  to  form  hyperspherical  clusters  is  a  major  limi¬ 
tation.  Due  to  the  temporal  dimension,  observations  associated  with  the  correct 
cluster  will  form  tracks,  not  hyperspheres.  This  result,  therefore,  necessitates  the 
development  of  new  clustering  methods  which  explicitly  account  for  the  temporal 
dimension  in  order  to  be  useful  in  the  multi-target  tracking  problem. 

1.5  Current  Approach 

Due  to  the  inherent  complexities  of  the  multi-target  tracking  problem, 
it  is  a  formidable  task  to  develop  an  algorithm  which  can  be  shown  to  be  optimal 
and  possess  the  following  characteristics, 

•  Perform  both  track  initiation  and  track  maintenance  and 

•  Permit  processing  of  data  in  real  time  while  minimizing 

-  Computational  complexity  and 

-  Data  storage  requirements. 

The  last  two  sub-objectives  are  important  not  only  in  achieving  real  time  per¬ 
formance  but  also  in  simplifying  the  processing  component  of  the  space-based 
sensor. 

A  heuristic  method  is  developed  which  combines  the  most  attractive 
features  of  the  non-hierarchical  clustering  approaches  with  the  track  initiation 
and  track  maintenance  approaches  suggested  in  references  [38,39]  (Morefield), 
[16]  (Chang  and  Youens),  [17]  (Chang  and  Tabaczynski),  and  [13]  (Blackman). 
In  fact,  Blackman  [13]  addresses  the  idea  of  combining  these  features  in  the 
track  maintenance  phase,  but  does  not  provide  a  workable  track  initiation  process 
capable  of  handling  large  numbers  of  targets  in  ballistic  trajectories.  Without 
track  initiation,  track  maintenance  cannot  be  performed. 

To  demonstrate  the  effectiveness  of  this  heuristic  approach,  the  method 
developed  is  specifically  tailored  to  the  ballistic  missile  defense  problem  with  tar¬ 
gets  (ICBMs)  in  flight  above  a  spherical  earth  with  no  atmosphere.  Orbiting 
satellite  sensors  surveying  the  ICBM  attack  provide  time-tagged  observations  of 


each  target’s  range,  range  rate,  azimuth,  and  elevation.  Additional  details  regard¬ 
ing  the  simulation  design  and  modeling  assumptions  are  provided  in  Chapter  4 
and  Appendix  A. 

While  the  track  initiation  process  marks  the  beginning  of  a  track’s  life 
cycle,  the  temporal  clustering  process  itself  begins  with  the  track  maintenance 
phase.  As  seen  in  Figure  1.2,  the  temporal  clustering  process  begins  by  reading 


Figure  1.2:  Temporal  Clustering  Process 

in  the  data  from  the  current  observation  frame  and  then  forecasting  all  existing 
clusters  to  the  current  observation  time.  The  costs  of  associating  each  new  ob¬ 
servation  with  a  predicted  observation  corresponding  to  each  existing  cluster  are 
calculated  and  observations  satisfying  the  gating  criteria  are  assigned  to  clusters 
to  minimize  the  total  overall  association  costs.  Each  cluster  receiving  a  new  ob- 
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servation  has  its  state  estimate  updated  while  clusters  not  receiving  observations 
are  considered  for  termination.  Finally,  all  remaining  observations  are  passed 
to  the  track  initiation  procedure  for  evaluation  in  determining  feasible  sets  of 
observations  to  initiate  new  tracks. 

In  Chapter  2,  the  various  components  making  up  the  track  maintenance 
process  of  assigning  observations  from  the  current  observation  frame  to  existing 
clusters  will  be  discussed.  The  process  of  deciding  how  and  when  to  terminate  a 
track  is  also  addressed  here. 

Then,  in  Chapter  3,  an  effective  means  for  performing  track  initiation 
is  developed.  Chapter  4  contains  a  discussion  of  the  specific  simulation  scenarios 
examined  as  well  as  an  analysis  of  the  results.  Finally,  conclusions  are  presented 
in  Chapter  5  together  with  a  discussion  of  proposed  extensions  to  the  current 
research. 


Chapter  2 


Track  Maintenance 


Each  observation  frame  read  by  the  sensor,  yields  a  set  of  observation 
vector^  consisting  of  the  observation  time  and  each  object’s  range,  range  rate, 
azimuth,  and  elevation  relative  to  the  sensor.  This  set  of  vectors  does  not  neces¬ 
sarily  contain  observations  of  all  the  targets  in  the  sensor’s  field-of-view.  This  is 
due  to  observability  problem  ising  from  sensor  characteristics  and/or  defects  or 
as  a  result  of  the  sensor 1  arget  geometry.  In  addition,  there  may  be  observations 
which  do  not  correspond  io  any  physical  target,  but  are  again  the  result  of  sensor 
characteristics  and  the  observation  environment.  Many  of  these  spurious  mea¬ 
surements  can  be  eliminated  by  pre-processing  the  data  to  remove  inconsistent 
observations  in  light  of  the  sensor  characteristics. 

Since  the  track  maintenance  process  is  restricted  to  dealing  with  only 
those  observations  recorded  by  the  sensor,  it  must  be  capable  of  ;  igning  ob¬ 
servations  to  existing  clusters  so  that  inappropriate  assignments  are  disallowed 
without  eliminating  correct  assignments.  This  capability  is  necessary  to  prevent 
making  assignments  to  a  track  when  observations  are  missing  from  that  track.  It 
must  also  be  able  to  continue  existing  tracks  which  do  not  receive  an  assignment 
until  such  time  as  track  termination  is  deemed  appropriate. 

2.1  Forecasting 

To  decide  whether  an  observation  should  be  considered  for  assignment 
to  an  existing  cluster,  the  “closeness”  of  each  observation  to  each  cluster  must 
first  be  determined.  To  do  so,  however,  both  the  observation  and  the  cluster  must 
be  evaluated  in  the  same  space  and  at  the  same  epoch.  In  the  problem  under 
consideration,  two  spaces  are  used:  the  attribute  space  and  the  state  space. 
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In  this  study,  the  attribute  space  is  a  four-dimensional  spherical  coor¬ 
dinate  space  centered  on  the  satellite  sensor.  The  direction  of  its  primary  axis 
is  fixed  and  points  toward  the  vernal  equinox.  The  four  dimensions  are  the  four 
attributes  measured  for  each  target:  range,  range  rate,  azimuth,  and  elevation. 

An  object’s  state  is  some  set  of  physical  characteristics  which,  together 
with  a  knowledge  of  the  state  transition  rules,  allows  predictions  of  the  state  at 
any  future  time.  For  this  study,  the  state  space  is  a  six-dimensional  inertial  Carte¬ 
sian  coordinate  space  centered  at  the  center  of  mass  of  the  Earth.  The  direction 
of  its  primary  axis  also  points  toward  the  vernal  equinox.  This  coordinate  system 
is  referred  to  as  the  Earth  Centered  Inertial  (ECI)  coordinate  system.  The  six 
dimensions  of  this  space  consist  of  three  rectangular  position  components  and 
three  rectangular  velocity  components. 

The  choice  of  this  state  space  is  possible  because  the  motion  of  a  target 
in  a  conservative  force  field  (one  which  can  be  described  completely  by  a  poten¬ 
tial)  can  be  fully  described  given  that  target’s  position  and  velocity.  The  case  un¬ 
der  examination  includes  only  gravitational  effects  and  excludes  non-conservative 
forces  such  as  thrust  and  drag,  and  is,  therefore,  a  conservative  system. 

The  use  of  the  state  space  has  several  advantages  over  that  of  the  at¬ 
tribute  space  in  the  temporal  clustering  process.  Data  storage  is  minimized  be¬ 
cause  knowledge  of  a  target  can  be  maintained  in  a  single  state  vector  rather  than 
a  track  of  observations.  Methods  for  efficiently  tracking  targets  in  the  state  space, 
such  as  the  Kalman  filter,  are  readily  available.  And,  for  this  study,  estimates 
in  the  state  space  can  be  easily  and  unambiguously  mapped  into  the  attribute 
space  while  the  converse  is  not  true. 

To  begin  the  process  of  assigning  observations  to  clusters,  therefore, 
each  cluster’s  state  vector  is  projected  to  some  common  observation  epoch  and 
then  mapped  into  the  attribute  space.  Not  only  is  it  necessary  to  forecast  the 
state  vector  to  the  observation  epoch,  however,  but  the  associated  state  covariance 
matrix  to  be  used  to  gate  the  observations  must  also  be  forecast  to  that  epoch  and 
both  the  state  and  state  covariance  matrix  must  be  mapped  into  the  attribute 
space  for  direct  comparison  with  the  observations.  As  seen  in  Figure  2.1,  the 
mapped  state  covariance  can  be  used  to  form  a  confidence  interval  around  the 
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projected  state  estimate.  Two  tracks  are  shown  together  with  their  estimated 
states  and  gates.  Only  observations  (hollow  circles)  falling  within  a  gate  are 
considered  for  possible  assignment  to  an  existing  track.  In  this  example,  one 
observation  could  be  assigned  to  either  Track  A  or  B,  one  observation  could  be 
assigned  only  to  Track  B,  and  one  observation  could  be  assigned  to  neither. 


Track  A 


Track  B 

Figure  2.1:  Gating  Process 

The  original  state  covariance  matrix  is  determined  when  the  initial  state  estimate 
is  formed  and  is  described  in  detail  in  Section  3.3.1. 

2.1.1  Kalman  Filter 

A  natural  method  for  forecasting  the  cluster  state  vector  is  provided  by 
the  Kalman  filter.  Not  only  does  the  Kalman  filter  provide  a  recursive  means  of 
propagating  the  state  estimate  and  state  covariance  matrix  but  it  also  provides  an 
optimal  means  for  updating  the  same  with  the  current  observation.  The  Kalman 
filter,  as  developed  in  Kalman  [33]  and  Kalman  and  Bucy  [34],  assumes  that 
discrete  states  are  linearly  related  via  a  state  transition  matrix  and  that  discrete 
observations  are  linearly  related  to  the  current  state.  That  is 

S*+i  =  $(<jt+i>  <jfe)Sfc  +  Ufc  (2.1) 

O*  =  HfcSfc-fw*,  (2.2) 


where 


£[u*]  =  0  SKu^]  =  Q  kSjk 
£[w*]  =  0  E{  WjWjH  =  R  k6jk 


£[ujW]  =  0. 
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Here  Sjt  is  the  target  state  at  time  t*,  tk)  is  the  state  transition  matrix 

which  governs  how  S  changes  from  time  tk  to  tfc+i,  and  u*  is  a  white  noise  process. 
O k  is  the  measurement  for  the  target  at  time  tk,  Hk  is  the  measurement  transform 
matrix  which  governs  how  the  state  Sk  and  the  measurement  O*  are  related,  and 
Wfc  is  another,  independent,  white  noise  process.  Both  Qt  and  Rfc  are  assumed 
diagonal  matrices  and  6jk  is  the  Kronecker  delta. 

The  optimality  criterion  is  that  the  estimate  be  a  minimum  variance 
unbiased  estimate  of  the  true  state.  That  is 


Minimize  E[e£  R^e*] 

subject  to  i?[Sfc]  =  Sjt, 

(2.4) 

(2.5) 

where 

e*,  =  (Ofc  -  HfcSjfe). 

(2.6) 

Assuming  prior  estimates  of  the  state  Sk  and  the  state  error  covariance 
Pjt,  the  resulting  minimum  variance  unbiased  estimate  is 

Si*  =  S/t  +  Kfc(Ofc  —  HfcSfc), 

(2.7) 

where 

Kfc  =  P*;Hf(HjtPfcHfc  +  Rfc)_1- 
The  updated  state  error  covariance  then  becomes 

(2.8) 

Pfc  =  (I  -  KfcHfc)P* 

(2.9) 

and  the  state  and  state  error  covariance  matrix  are  propagated  according  to 

Sfc+1  =  *(**+  (2.10) 

Pfc+l  =  *(tk+l,tk)Pk$(tk+1,tk)T  +  Qk.  (2.11) 

Unfortunately,  for  a  target  in  a  ballistic  trajectory,  neither  the  dynamics 
nor  the  measurements  are  linear.  However,  the  system  can  be  linearized  through 
the  application  of  Taylor  series  expansions. 


2.1.2  Extended  Kalman  Filter 


For  a  system  with  nonlinear  dynamics  and  nonlinear  measurements,  the 
standard  model  of  Equations  2.1  and  2.2  becomes 

S  =  f(S,t)  +  u(<)  (2.12) 

O  =  h(S,<)  +  w(<)  (2.13) 

with  the  covariance  properties  of  u(f)  and  w(f)  unchanged.  Equation  2.12  is  the 
differential  equation  describing  how  the  system  dynamics  affect  the  state,  S,  and 
Equation  2.13  describes  the  specific  relationship  between  the  observation  and  the 
state  over  time. 

Given  some  nominal  reference  trajectory  S *(t),  the  true  trajectory  may 
be  written  as 

S(0  =  S*(i)  +  s(0  (2.14) 

so  that  Equations  2.12  and  2.13  become 

S*  +  s  =  f(S*  +  s,  f )  +  u(t)  (2.15) 

O  =  h(S*  +  s,f)  +  w(f).  (2.16) 

Assuming  s  to  be  small,  f  and  h  may  be  approximated  with  Taylor 

series  expansions,  so 

S*  +  s  =  f(SV)  +  (|0  S  T  •  •  •  +  u(0  (2.17) 

O  =  h(S*,f)+ sT-'-Tw (0,  (2.18) 

where 

S*  =  f(S*,f),  (2.19) 


Retaining  only  first-order  terms,  a  linearized  model  in  s  results  as 

S  =  S  +  U(<)  (2.21) 

°  =  (^|)  S  + w(t)  =  ° -h(s*,0  (2.22) 

or,  defining 

A«)=(|)'  and  H  =(|)‘,  (2.23) 

then 

s  =  A(f)s  +  u  (t)  (2.24) 

o  =  Hs  +  w(f).  (2.25) 

Standard  algorithms  for  implementing  the  Kalman  filter  can  be  used, 
with  the  only  difference  being  that  the  nominal  state  vector,  S’,  and  the  state 
covariance  matrix,  $(4+1,4),  are  updated  through  the  use  of  a  numerical  inte¬ 
gration  routine  between  time  steps,  with 

S*  =  f(SV)  $(*,  4-i)  =  A(0$(44_i) 

and  (2.26) 

S*(4-i)  =  SJ_,  $(4-i,  4-i)  =  I- 

This  implementation  is  known  as  the  linearized  Kalman  filter. 

Frequently  it  is  desirable  to  use  the  current  estimated  trajectory  in 
place  of  the  nominal  reference  trajectory  since,  under  stable  conditions,  a  better 
estimate  will  result.  The  process  of  updating  the  reference  trajectory  with  the 
latest  estimate  of  the  true  trajectory  is  known  as  rectification.  Using  rectification, 
the  reference  trajectory  is  updated  as 

SI  =  SI  +  sk  (2.27) 


for  each  new  measurement.  This  rectified  linearized  Kalman  filter  is  more  com¬ 
monly  known  as  the  Extended  Kalman  Filter  (EKF)  and  is  the  filter  of  choice 
for  most  astrodynamical  tracking  problems. 


2.1.3  State  Propagation 

At  time  f*.,  the  state  vector,  S^,  consists  of  the  target’s  position,  rjt,  and 
velocity,  f *  (or  v*),  in  ECI  coordinates 


Sit  =  {xk,yk,Zk,ik,yk,Zk)T- 


The  observation  vector,  Ojt,  is 


O/e  —  (/^fci  pk:  &k)  i 


(2.28) 


(2.29) 


where  pk  is  the  range  from  the  sensor  to  the  target,  pk  is  the  range  rate,  or*  is 
the  target’s  azimuth,  and  £*  is  the  target’s  elevation. 

Given  estimates  for  Sm  and  Pm  at  some  time  tm  prior  to  the  time  of 
the  current  observation,  the  method  to  be  used  to  forecast  the  target  state  and 
covariance  to  the  current  observation  time,  t0,  can  be  developed  using  the  specific 
dynamical  model  for  this  investigation. 

As  shown  in  Equation  2.26,  the  target  state  estimate,  S,  is  updated 
through  numerical  integration  of 


with  initial  conditions 


S  =  f(S,f) 


S  (tm)  =  S„ 


(2.30) 


(2.31) 


to  the  current  observation  time,  tQ.  The  specific  function  f(S,f)  depends  on  the 
system  dynamics.  In  general,  for  a  target  in  ballistic  flight, 


S  =  G(r,t)  +  D(r,  v,f)  +  +  71(f). 


(2.32) 


Each  term  on  the  right-hand-side  of  Equation  2.32  represents  a  force  per  unit 
mass  (i.e.,  acceleration)  on  the  target.  G(r,  t)  is  the  combined  gravitational 
acceleration  on  the  target  due  to  the  earth,  sun,  and  moon,  D(r,  v,  t )  is  the  effect 
of  atmospheric  drag,  T(f)  is  the  target’s  thrust,  M(t)  is  the  target’s  mass,  and 
7Z(t)  is  the  acceleration  due  to  all  unmodeled  forces.  For  this  study,  the  target  is 
subjected  only  to  the  gravitational  acceleration  of  a  spherical  earth. 
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For  such  a  target  in  a  two-body  orbit,  the  standard  second-order  differ¬ 
ential  equation  is 

S  =  (2.33) 

which,  when  expressed  in  the  first-order  form  of  Equation  2.12  gives 


S  =  f(S,t)  =  ( x,y,z 


r3  ’  r3  , 


(2.34) 


and  can  be  integrated,  using  the  initial  condition  S (tm)  =  Sm,  from  tm  to  tQ. 
Equation  2.34  results  from  differentiating  Equation  2.28  and  applying  Equa¬ 
tion  2.33. 

Propagating  the  state  covariance,  Pm,  is  a  bit  more  difficult.  From 
Equation  2.11, 

PD  =  *(t0,tm)Pm*(t0,tm)T  +  Q,  (2.35) 

where  the  white  noise  process  u  is  assumed  to  be  homoscedastic  and  known. 
Otherwise  adaptive  filtering  techniques  are  necessary  to  adaptively  compute  Qm. 
But,  4 >(t0,tm)  itself  must  be  numerically  integrated,  according  to  Equation  2.26 
and  the  specific  form  of  A(t)  must  be  evaluated. 

Letting  S  =  (r,v)T  and  applying  the  definition  of  A (t)  from  Equa¬ 
tion  2.23, 

«»- (8)’ -(!)-(': £)■ 


where 


.  dr 

frr  “  Fr  ~  ’ 

f  =—  =  1 

rv  dv  ’ 


f  _  Hi  - 


(2.37) 


(2.38) 


‘-3(f)2  -3(ff)  -3(f)  \ 

-3(f)  1-3(8)’  -3(f)  .(2.39, 

-3(ff)  -3(f)  .-3(f)2] 


and  fvv  =  7T-  =  0. 
9v 


(2.40) 
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Therefore,  the  system  of  equations  to  be  integrated  is 


d>  $ 

^rr  ^rv 


subject  to 


So,  to  find  PQ,  3>(t0,tm)  is  integrated  according  to  Equations  2.42  through  2.47 
and  Equation  2.35  is  applied. 

With  S0  and  PQ  the  state  estimate  and  its  covariance  must  now  be 
mapped  into  the  attribute  space  for  use  in  completing  the  cluster  assignment 
process. 

2.1.4  Mapping  Into  the  Attribute  Space 

Given  the  target  state,  S0  =  (r0,r0)T,  and  the  sensor  state,  Sos  = 
(roa?ro»)Ti  the  transformation  from  the  state  space  into  the  attribute  space  is 
given  by 


Po  (Poxi  Poyi  Poz)  —  rOJ 

(2.48) 

Po  ( Pox  i  Poy  >  poz )  =  —  ros 

(2.49) 

po  =  IIPoll 

(2.50) 

po  =  Po  Po 

(2.51) 

a0  =  tan-1  (— ) 

\Pox) 

(2.52) 

If 

I 

H 

I 


0  l\/$„  *rv  \ 

Ur  0  )  \  $ur  / 

(2.41) 

■  K|'i 

=  S.r 

=  &vv 

(2.42) 

(2.43) 

1 

-  f  * 

^rr 

=  fvr&rv 

(2.44) 

(2.45) 

r 

)  ^m(^mi  tm)  I 

(2.46) 

)  =  ^tir(tm>tm)  =  0 

(2.47) 

1 

'» ■.*** ,*<•'** t  '«»v  i 


_rfB  ^tb  *  *  #<>  »-.,-  t - 1 


■  l*i  I  jj**  i 


I  .%'  •■*•'  »_«>*¥«,*•  .**§  •  (•-*'  -^1  ? .* 


e0  =  tan 


-i 

VftJ  ’ 


where 


0o  —  V  Pox2  +  Poy2 


(2.53) 


(2.54) 


and  p0  is  the  unit  vector  along  the  range  vector,  pQ.  These  transformations  are 
based  upon  the  sensor-target  geometry  depicted  in  Figure  2.2  and  a  standard 
transformation  of  rectangular  coordinates  of  the  state  space  into  the  spherical 
coordinate  system  of  the  attribute  space. 


Target 


Sensor 


Figure  2.2:  Sensor-Target  Geometry 

The  vector  of  estimated  observations  at  the  current  observation  time, 
<0,  is  designated  by  00,  where 


G0  (Poi  Pot  &o)  G^Sq). 


(2.55) 


More  properly,  00  should  be  considered  to  be  a  function  of  the  uncertain  param¬ 
eters  r0  and  r0,  with  ros,  roa,  and  t0  treated  as  constants. 

Using  a  Taylor  series  expansion  and  retaining  only  first-order  terms,  the 
state  estimate  can  be  shown  to  be  an  unbiased  estimate.  That  is 


£[o0]  =  e[60(  s0)] 


=  E  °o(»S J  + 


(S0  -  Pso) 


(2.56) 

(2.57) 
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where 


=  W  +  (gf)  *S.]-g£)  «. 

=  O0(^so)> 


(dOoV  d60 

UsJ  dS0s0=M(l 


(2.59) 

(2.60) 


(2.61) 


Now,  to  find  the  state  covariance  matrix,  a  Taylor  series  expansion  is  again  ap¬ 
plied,  so 


o0  -  E[o0]  =  60(fi$J  +  H0(S0  -  usj  +  •  •  •  -  E[6a] 
=  H0(S0  —  /i§o)  +  •  •  ■ 


and,  retaining  only  first-order  terms, 


s .  =  b  |(6,  -  E{ 60])(6„  -  f;[6j)T] 

=  E  [(H0(So  -  ^s.))(H„(S„  -  msJ)T] 
=  E  [(H„(S0  -  „sJ(S„  -  MS.)THr] 

=  H0£  j(S„  -  MS.MS.  -  xs.)Tj 

E  [(§„  -  /.S„)(S„  -  Ms.)r] 


(2.62) 

(2.63) 


(2.64) 

(2.65) 

(2.66) 
(2.67) 


(2.68) 


is  merely  the  state  covariance  matrix,  Pc.  So,  the  estimated  attribute  covariance 
matrix,  30,  can  be  written 

30  =  H0P0H„ .  (2.69) 

Actually,  30  =  max(H0P0H^,  R)  is  used,  where  R  is  the  observation  covariance 
matrix.  This  prevents  the  estimated  error  associated  with  the  observation  from 
becoming  smaller  than  the  known  error  of  the  sensor. 


I  i>i  (i.  it.  I'l'iVlK'iWiL' 


Also,  since  /igo  is  not  known,  the  reference  state,  S*,  is  used  instead,  so 


nQ  = 


dS0  lso=S;' 


(2.70) 


In  particular, 


^  0  0  0 

Po 


H0  = 


PoxPo  PoxPo  Poypo  PoyPo  PzPo  PozPo  Pox  Poy  Poz 
p?  P2  Pc2  Po  Po  Po 


(2.71) 
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Once  reasonable  estimates  of  Oc  and  30  have  been  determined,  they 
can  be  used  to  gate  the  new  observations. 

2.2  Cluster  Assignment 

Now,  each  actual  observation  is  compared  to  each  predicted  observation 
and  the  “cost”  of  association  (i.e.,  the  cost  of  assigning  an  actual  observation  to 
an  existing  cluster)  is  computed.  This  cost  is  based  upon  the  Euclidean  metric 
in  the  attribute  space.  Because  of  the  disparate  scales  of  the  various  attributes, 
each  attribute  of  the  observations  is  standardized  by  subtracting  the  attribute 
minimum  and  dividing  by  the  attribute  range. 

To  minimize  the  likelihood  of  infeasible  associations,  actual  observations 
which  fall  outside  the  predicted  observations’  gates  are  assigned  an  arbitrarily 
large  cost.  These  gates  are  derived  from  the  diagonal  elements  of  the  attribute 
covariance  matrix,  a0.  While  the  use  of  simple  rectangular  gates  based  on  the 
diagonal  covariance  elements  may  be  conservative  (depending  upon  the  relative 
magnitude  of  the  off-diagonal  terms),  it  is  shown  empirically  to  work  quite  well. 
An  association  is  considered  to  be  infeasible  if  any  single  attribute  of  an  ob¬ 
servation  is  outside  the  predicted  attribute’s  6-cr  confidence  interval  or  if  any 
two  attributes  are  outside  the  corresponding  predicted  attributes’  3-cr  confidence 
intervals. 


The  rationale  for  using  a  two-step  voting  process  for  determining  wheth¬ 
er  an  association  is  infeasible  or  not  is  based  upon  the  probability  of  Type  I  and 
Type  II  errors.  In  all  the  simulation  scenarios  run  to  date,  each  incorrect  associ¬ 
ation  resulted  in  at  least  one  attribute  being  far  outside  the  predicted  attribute’s 
3-<t  confidence  interval,  therefore  making  the  probability  of  a  false  association 
(Type  II  error),  even  at  the  6-<r  confidence  level,  quite  small. 

And  while  the  probability  of  denying  a  correct  association  (Type  I  error) 
at  the  3-<t  confidence  level  is  small,  the  likelihood  of  such  an  occurrence  over  the 
life  of  the  scenario  must  be  considered.  The  probability  of  a  Type  I  error  is 
equal  to  the  probability  that  one  or  more  observation  attributes  fall  outside  their 
gates,  or  one  minus  the  probability  that  none  fall  outside  their  gates.  If  p  is  the 
probability  that  a  single  observation  attribute  is  within  its  3-cr  gate  and  q  =  1  —  p, 
and  n  is  the  number  of  attributes,  then  the  probability  of  disallowing  a  correct 
assignment  is 


P(gating  error)  =  1 


(;) 


pV 


i-pn. 


(2.72) 


Applying  Equation  2.72  to  a  typical  100- target  scenario  with  four  attribute  types, 
1.08  Type  I  errors  would  be  expected  in  each  observation  frame,  resulting  in  a 
lost  observation.  By  requiring  that  two  attributes  (out  of  four)  be  outside  the 
3-cr  confidence  interval  the  probability  of  a  Type  I  error  is  reduced  considerably 
(by  a  factor  of  250). 

Once  the  costs  of  association  have  been  computed,  an  assignment  can 
be  made  of  actual  observations  to  existing  clusters  so  that  the  total  cost  of  these 
associations  is  a  minimum.  First,  obvious  assignments  are  made  where  only 
one  assignment  is  possible.  Then,  since  the  remaining  numbers  of  clusters  and 
observations  will  likely  be  unequal,  dummy  clusters  or  observations  are  formed 
with  association  costs  set  to  an  arbitri.  ,!y  large  value.  The  remaining  assignment 
is  then  solved  (in  this  case  using  the  Hungarian  algorithm  [40]),  and  all  clusters 
receiving  legitimate  assignments  are  updated  via  the  EKF  as  shown  in  Section  2.4. 


All  clusters  for  which  no  observation  can  feasibly  be  assigned  are  consid¬ 
ered  to  have  missed  an  observation  and  are  annotated  to  indicate  this.  Obviously, 
no  updating  of  this  cluster’s  state  or  state  covariance  is  possible.  Finally,  all  re¬ 
maining  unassigned  observations  are  then  passed  to  the  track  initiation  algorithm 
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discussed  in  Chapter  3. 

2.3  Track  Termination 

Conditions  for  terminating  a  cluster  are  now  considered.  Target  tracks 
may  require  termination  because  they  have  reached  their  impact  point,  been  de¬ 
stroyed,  become  obscured  by  the  earth’s  surface  or  atmosphere,  or  simply  because 
they  exit  the  sensor’s  field-of-view. 

In  the  scenarios  investigated  in  this  study,  all  observation  frames  are 
recorded  at  equally  spaced  intervals  and  all  observations  in  the  same  frame  have 
the  same  time  tag,  it  makes  sense  to  terminate  a  cluster  after  some  pre-specified 
number  of  missing  observations.  There  are  several  practical  reasons  for  so  doing. 
First,  there  is  no  point  (computationally)  in  continuing  to  propagate  a  cluster 
which  has  either  disappeared  or  been  terminated  due  to  a  Type  I  error.  After 
some  period  of  time  it  must  be  accepted  that  the  target  is  no  longer  being  tracked. 
With  fixed  time  interval  observation  frames  a  maximum  number  of  consecutive 
missing  observations  can  be  used  as  a  limit. 

Another  reason  for  not  propagating  clusters  indefinitely  relates  to  the 
increased  probability  of  Type  II  errors.  As  a  cluster  is  propagated  without  up¬ 
dating  its  state  and  state  covariance,  the  state  covariances  will  grow,  making  it 
more  and  more  likely  that  a  false  association  will  result.  This  limitation  could 
lead  to  another  criterion  for  terminating  a  cluster.  That  is,  the  cluster  could  be 
terminated  when  the  state  covariance  elements  exceed  certain  bounds.  The  main 
drawback  to  this  criterion,  however,  is  that  the  magnitude  of  the  covariance  el¬ 
ements  are  highly  dependent  upon  the  sensor-target  geometry,  so  determination 
of  the  bounds  would  be  subjective. 

In  using  the  consecutive  missing  observations  limit  as  a  criterion  for 
terminating  a  cluster,  the  likelihood  of  incorrectly  terminating  an  active  cluster 
due  to  a  chance  occurrence  of  the  limit  being  exceeded  must  be  considered.  To 
do  this,  an  upper  bound  on  the  likelihood  of  an  observation  being  missing  must 
be  estimated.  If  the  maximum  number  of  missing  observations  allowed  is  n  and 
the  probability  of  an  observation  being  missing  is  q  (assuming  equally  likely  and 
independent  events),  and  p  =  1  —q,  then  the  probability  of  incorrectly  terminating 
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a  cluster  is 


P(false  termination)  =  +  ^p°gn+1  =  gn+1 .  (2.73) 


For  a  100-target  scenario  with  5-second  data  intervals  covering  the  first 
six  minutes  of  an  ICBM  launch,  the  following  number  of  false  terminations  based 
upon  chance  are  expected: 

Maximum  Probability  of  Missing 


Missing 

0.05 

0.10 

0.20 

1 

365.00 

730.00 

1460.00 

2 

18.00 

72.00 

288.00 

3 

0.89 

7.10 

11.20 

4 

0.04 

0.70 

2.21 

5 

0.00 

0.07 

0.44 

6 

0.00 

0.01 

0.09 

7 

0.00 

0.00 

0.02 

Table  2.1:  Expected  Number  of  False  Terminations 

Other  criteria  for  terminating  clusters  could  consider  whether  the  target 
was  predicted  to  be  beyond  the  field-of-view  of  the  sensor,  below  the  surface  of 
the  earth,  or  beyond  the  earth’s  horizon.  However,  these  considerations  are  not 
presently  implemented  in  the  simulation. 

2.4  Track  Update 

Once  an  assignment  has  been  made  between  the  existing  clusters  and 
the  observations  at  the  current  observation  time,  tQ,  the  state  for  each  cluster  is 
ready  to  be  updated.  Using  the  linearized  versions  of  Equations  2.7,  2.8,  and  2.9, 


where 


s0  =  s0  +  K0(oo  —  H0s0) 
and  P0  =  (I  -  K0H0)Po, 

K0  =  P0H^(HoP0H^  +  R)_1. 


(2.74) 

(2.75) 

(2.76) 


n 

SC 


Again,  R  is  used  rather  than  R0  in  Equation  2.76  under  an  assumption  that  the 
white  noise  process  w(f)  is  homoscedastic.  And,  because  the  Extended  Kalman 
Filter  is  being  used,  s0  =  0,  thus,  Equation  2.74  becomes  (assuming  rectification 
occurs  at  each  state  update) 

So  =  KoOo,  (2.77) 


where 


o0  =  o0  -  o;. 


(2.78) 


The  final  item  needed  is  the  measurement  transform  matrix,  H0  from 
Equation  2.71.  Evaluating  H0  at  S0  and  forming  K0  according  to  Equation  2.76, 


S0  —  Ro^o) 

So  =  So  +  So, 

and  P0  =  (I  -  KoHo)P0. 


(2.79) 

(2.80) 
(2.81) 


While  the  implementation  discussed  above  is  theoretically  correct,  mod¬ 
ification  of  the  EKF  is  often  necessary  to  prevent  filter  divergence.  One  of  the 
primary  causes  of  filter  divergence  is  associated  with  errors  which  occur  in  the 
computation  of  the  state  error  covariance  matrix  [55,25].  In  particular,  round¬ 
off  errors  in  the  calculation  of  this  matrix  can  cause  it  to  become  non-positive 
semi-definite — a  theoretical  impossibility. 

In  this  investigation,  the  standard  EKF  formulation  was  found  to  suffer 
from  just  this  type  of  filter  divergence.  The  state  error  covariance  matrix  immedi¬ 
ately  became  non-symmetric  and  non-positive  semi-definite  during  the  first  filter 
update.  This  result  was  due  to  the  extremely  poor  conditioning  of  the  matrix 


M  =  HPHr  +  R, 


(2.82) 


which  must  be  inverted  in  Equation  2.76  as  part  of  the  EKF'  procedure.  As  a 
result,  the  EKF  is  re-formulated  to  take  advantage  of  an  approach  which  main¬ 
tains  the  natural  symmetry  and  positive  semi-definiteness  of  the  state  covariance 
matrix. 


2.4.1  Matrix  Square  Root 

Since  the  a  priori  state  error  covariance  matrix  P  is  a  symmetric  positive 
semi-definite  matrix,  it  can  be  written  as 

P  =  WWr,  (2.83) 


where  W  is  the  matrix  square  root  of  P.  Such  a  square  root  can  easily  be  found 
using  Cholesky’s  Decomposition  Algorithm  [55].  Given  the  n  x  n  elements  of  P, 
the  elements  of  W  may  be  found  using  the  following  procedure 


For  i  =  1,2, ...  ,n 


p«  -  e  n 

k=  1 


For  j  =  *  +  !,...,» 


(2.84) 


Pij  - 

^ - •  (2-85) 

The  resulting  matrix  is  a  lower-triangular  square  root  matrix.  The  state  error 
covariance  update  can  now  be  reformulated  using  W  instead  of  P  to  ensure  that 
P  remains  symmetric  and  positive  semi-definite. 


2.4.2  Covariance  Update  Reformulation 

From  Equations  2.75  and  2.76,  the  state  error  covariance  update  equa- 


tion  is 

P  =  (I  -  KH)P, 

(2.86) 

where 

K  =  PHr(HPHT  +  R)-\ 

(2.87) 

Equivalently, 

P  =  P  -  PHt(HPHt  +  R)_1HP. 

(2.88) 

Substituting  WWT  for  P  and  WWT  for  P  yields 
WWT  =  WWT  -  WWtHt(HWWtHt  +  R)1HWWt. 


(2.89) 


Letting  F  =  WTHT  and  M  =  FTF  +  R,  then  Equation  2.89  becomes 


WWT  =  W(I  -  FM“1Ft)Wt. 

(2.90) 

Expressing  (I  —  FM-1Fr)  as  AA7,  then 

ww7  =  WAA7W7 

(2.91) 

or 

W  =  WA. 

(2.92) 

The  state  error  covariance  update  is  found  by  first  computing 

* 

II 

*0' 

(2.93) 

using  Cholesky’s  Decomposition  Algorithm,  then  forming 

F  =  WtHt 

(2.94) 

M  =  FtF  +  R 

(2.95) 

A  =  (I  —  FM_1Ft)2 

(2.96) 

and  W  =  WA 

(2.97) 

Since  FTF  -f  R  is  symmetric,  symmetric  inverse  routines  can  be  used  to  calculate 

M"1. 

Obviously,  P  =  WWT  and 

s  =  s  +  K(o  -  Hs), 

(2.98) 

where 

K  =  WFM"1, 

(2.99) 

which  for  the  EKF  becomes 

s  =  Ko. 

(2.100) 

This  method  ensures  that  P  remains  symmetric  and  positive  semi-defi¬ 
nite  as  expected. 
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Chapter  3 


Track  Initiation 


Once  the  track  maintenance  process  has  been  completed,  the  remaining 
unassigned  observations  are  passed  to  the  track  initiation  algorithm.  These  unas¬ 
signed  observations  are  most  likely  the  result  of  new  targets  which  may  appear 
in  the  sensor  field-of-view  because  they  were  just  launched,  were  just  deployed  as 
a  multiple  independently-targeted  reentry  vehicle  (MIRV),  entered  the  sensor’s 
field-of-view,  or  emerged  from  being  obscured  by  the  earth’s  surface  or  atmo¬ 
sphere. 

The  goal  of  this  process  is  to  form  an  initial  estimate  of  a  cluster’s  state 
and  state  covariance  using  the  minimum  number  of  observations.  Examination 
of  various  orbit  determination  techniques  in  [12,24,29,30]  has  shown  Laplace’s 
method,  using  three  observations  of  a  target’s  range,  range  rate,  azimuth,  and 
elevation,  to  be  the  most  appropriate  method  for  determining  an  initial  state  es¬ 
timate  in  this  study.  How  this  method  is  applied  will  be  shown  in  Section  3.3.1. 
However,  since  three  consecutive  observations  cannot  be  guaranteed,  the  unas¬ 
signed  data  must  be  stored  in  a  buffer  to  permit  forming  the  necessary  combina¬ 
tions  of  observations. 

As  in  the  track  termination  process,  the  likelihood  of  missing  an  ob¬ 
servation  will  determine  the  size  of  the  buffer  required.  If  a  buffer  covering  the 
last  m  observation  frames  is  used,  the  probability  of  failing  to  correctly  initiate  a 
cluster  will  be  the  probability  that  at  most  one  observation  of  the  target  associ¬ 
ated  with  that  cluster  occurs  in  the  last  m  —  1  observation  frames  given  that  an 
observation  has  been  detected  in  the  first  observation  frame.  That  is, 


^(initiation  failure)  = 


m  —  1 
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where  p  is  the  probability  of  detecting  an  observation  and  q  =  1  —  p. 

As  with  the  track  termination  process,  the  expected  number  of  failures 
to  initiate  a  cluster  for  various  buffer  sizes  and  probabilities  of  missing  data  is 
summarized  below  in  Table  3.1.  The  table  assumes  a  100-target  scenario. 


Buffer 

Probability  of 

Missing 

Size 

0.05 

0.10 

0.20 

3 

9.75 

10.00 

36.00 

4 

0.73 

2.80 

10.41 

5 

0.05 

0.37 

2.72 

6 

0.00 

0.05 

0.67 

7 

0.00 

0.01 

0.16 

Table  3.1:  Expected  Number  of  Track  Initiation  Failures 

Once  the  size  of  the  track  initiation  buffer  is  determined,  the  track 
initiation  process  can  be  analyzed.  Since  the  objective  in  this  process  is  to  form 
observation  triples  which  can  be  assessed  for  suitability,  it  is  desirable  to  form  all 
feasible  triples  and  select  among  these  for  the  “best”  overall  assignment. 

3.1  Problem  Formulation 

To  find  the  “best”  overall  assignment  for  the  track  initiation  problem 
requires  that  the  problem  to  be  solved  be  defined  specifically  as  well  as  in  what 
sense  the  solution  is  best.  Briefly,  the  problem  is  to  form  triples  of  observations 
(tracks)  such  that  no  observation  is  included  in  more  than  one  track  and  that  the 
system  dynamics  are  not  violated.  But  there  must  be  some  means  for  assessing 
the  “cost”  of  associating  an  observation  with  a  track.  Then,  the  problem  becomes 
to  choose  a  set  of  tracks  which  minimize  this  association  cost  subject  to  the 
restriction  that  no  observation  be  used  in  more  than  one  track  and  that  no  system 
dynamics  be  violated. 

A  measure  of  the  quality  of  the  tracks  of  targets  in  ballistic  trajectories 
is  a  function  of  the  total  specific  energy,  £.  High  quality  is  associated  with  low 
values  of  £.  The  designers  of  an  ICBM  booster  are  faced  with  providing  a  system 
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for  delivering  nuclear  warheads  to  their  targets  for  the  minimum  cost.  Since 
the  cost  is  energy,  it  is  reasonable  to  assume  that  each  warhead  will  follow  a  low 
energy  trajectory  to  its  designated  impact  point.  Combining  this  assumption  with 
the  fact  that  a  ballistic  trajectory  is  a  minimum  energy  path  in  the  gravitational 
potential  results  in  the  conclusion  that  any  misgrouping  of  observations  into  a 
track  will  require  a  higher  energy  orbit.  Therefore,  if  one  pair  of  observations 
defines  an  orbit  and  a  second  pair  defines  a  similar  orbit  (with  the  mid-point 
observation  common  to  each  pair)  there  should  be  no  change  in  energy  between 
the  two  orbits  if  the  three  observations  are  properly  grouped  (in  a  conservative 
force  field).  Since  imperfect  measurements  are  involved,  the  assumption  is  that 
this  change  is  smallest  for  properly  grouped  observations  than  otherwise. 

The  resulting  problem  can  be  formulated  as  a  binary  quadratic  program. 
Consider  m  observation  frames  and  np  unassigned  observations  in  each  frame 
p  =  1,2,  ...,m.  Observation  i  in  frame  p  is  linked  with  observation  j  in  frame 
q  if  and  only  if  the  variable  xipjv  =  1,  otherwise,  xtpjq  =  0.  The  cost  associated 
with  a  given  observation  triple  (pi,qj,rk)  (i.e.,  the  tth  observation  from  frame  p, 
the  jth  observation  from  frame  q,  and  the  kth  observation  from  frame  r),  ctpjqkr , 
is  equal  to  the  specific  energy  of  the  orbit  defined  by  the  observation  triple  if  the 
system  dynamics  are  satisfied,  and  equal  to  infinity  otherwise. 

Problem,  BQP  ( Binary  Quadratic  Program) 


Minimize 


subject  to 


m— 2  m— 1  m  nq  „r 

^2/C'Vi<ikTXii>jqXiqkr 
p=  1  7>P  r>q  i=l  j=  1  k= 1 

np  nr 

y  '.xiQkT  P  =  1, 2, . . .  ,m  2 

,=1  k=l  q  —  p  +  1 , . . . ,  m  —  1 

r  =  9+l,...,m 

j  —  1,2,...,  Tlq 


p  =  1,2, ...  ,m  —  2 
q  =  p-f  1, . . . ,  m  -  1 
J  —  1 , 2, .  • . ,  rtq 

p  =  1, 2, . . .  ,m  —  2 
q  =  p  +  1,. . .  ,m  —  1 

l  —  1,2,...,  Tip 


j= 1 
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n. 

XXi9*r  <  1, 

9  =  P+  1, 

...  ,m  —  1 

(3.6) 

}=1 

r  =  q+  1,. 

..,m 

k  =  1,2,.. 

.,nr 

XiPjq,XjqkT  binary, 

P  —  1,2,.. 

.  ,m  —  2 

(3.7) 

q  =  p  + 1,- 

...  ,m  —  1 

H 

+ 

II 

..,m 

i  =  l,2,... 

,np 

j  =  1,2,.. 

■,nq 

*  =  1,2,.. 

Equation  3.3  (conservation  of  flow)  ensures  that  if  observation  q3  is  paired  with 
some  observation  p;,  then  it  is  also  paired  with  some  observation  r*.,  thus  guar¬ 
anteeing  that  a  triple  is  always  formed.  And  Equations  3.3  through  3.7  ensure 
that  each  observation  is  assigned  to  at  most  one  track,  and  vice  versa. 

Since  integer  programs  can  be  difficult  to  solve  efficiently  and  quadratic 
integer  programs  even  more  so,  the  problem  is  reformulated  into  a  binary  linear 
program  by  setting 

l lipjqkr  =  XipjqXjqkr-  (3.8) 


This  variable  ensures  that  only  observation  triples  are  formed.  That  is,  obser¬ 
vation  i  in  frame  p,  observation  j  in  frame  q,  and  observation  k  in  frame  r  are 
linked  if  and  only  if  the  variable  yipjqkr  =  1,  otherwise,  yiPjqkr  =  0. 

Problem  BLP  (Binary  Linear  Program) 


Minimize 


subject  to 


m-2m-l  m  np  nq  nr 

X!  X  X  X/  X/  ^ \^CipjqkrVipjqkT 
p=l  q>p  r>q  t=l  j=\  k=  1 


(3.9) 


np 


^JJipjqkr  <  1, 

P  =  1,2,.. 

.  ,m  —  2 

i  =  l 

q  =  p  + 1,- 

..,m- 

r  =  9+  1,. 

..,m 

i  =  1,2,.. 

.,ng 

n. 

k  =  1,2,.. 

■  ,nr 

^  '.Viviakr  —  1 » 

p=  1,2,.. 

. ,  m  —  2 

i=i 

9  =  P+  1,  • 

...,771  — 

r  =  q+  i,- 

.  .  ,  777 

i=  1,2,... 

,  77p 

*  =  1,2,.. 

.  ,77r 

(3.11) 
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^Vipjqkr  <  1, 
k= 1 


Vipjqkr  binary, 


Lemma  1  Formulations  BQP  and  BLP  are  equivalent. 


p  =  1,2, . . . ,  m  —  2 

(3.12) 
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q  =  p  +  1, . ..  ,m  -  1 

r  =  q+  1, .. .  ,m 

2  — ■  1,2,...,  Tip 

£1 

j  =  1,2, ...  ,n. 

p  =  1,2, . . .  ,m  —  2 

(3.13) 

q  =  p+  1, .  ..,m  -  1 

r  =  q  +  1, .. .  ,m 

2  1,2,...,  72p 

4  1  ,  2,  .  .  .  ,  Tlq 

k  =  1,2, . . .  ,nr. 

■f 

PROOF:  The  objective  functions  of  both  programs  are  equivalent  by  the  definition 
of  Vipjqkr  ■  To  complete  the  proof,  any  solution  to  either  formulation  must  be  shown 
to  satisfy  the  constraints  of  the  other. 

All  feasible  solutions  to  BQP  satisfy  BLP: 

Multiplying  both  sides  of  Equation  3.4  by  Xjqkr  and  applying  Equation 
3.8  proves  that  Equation  3.10  is  satisfied.  From  Equation  3.3, 


np  nr 


X/x‘Pi? 

P 

=  1,2,.. 

.  ,m  —  2 

(3.14) 

t=i 

fc=i 

V 

=  P  +  1, 

. . .  ,m  —  1 

r 

=  9  +  1, 

. . . ,  m 

j 

=  1,2,.. 

■,ng, 

7lr 

which  implies  y^x,afcr 

<  1, 

9 

=  P+  1, 

...  ,m  —  1 

(3.15) 

fc=i 

r 

=  9  +  1, 

. . .  ,m 

j 

=  1,2,.. 

Again,  multiplying  Doth  sides  of  Equation  3.15 

by  %ipjq  and 

applying  Equation 

3.8  proves  that  Equation  3.12  is  satisfied.  Finally,  Equation  3.11  is  shown  to  be 
satisfied  by  noting  that 

n, 

^  1,  P=  1,2,  ...,m  -2  (3.16) 

J=1  q  =  p  +  1,  •  •  •  ,m  —  1 

i  =  1 , 2, . . . ,  np 
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Tlq 


&m2 

P=  1,2,.. 

.  ,m  —  2 

(3.17) 

J= 1 

9  =  P  +  1, 

. . . ,  m  —  1 

i  =  1,2,... 

■,np 

n. 

n. 

5 ~2x'pjq  xiqkr 

—  1 

P  —  1,2,.. 

. ,  m  —  2 

(3.18) 

j- 1 

q  =  p  + 1, 

. . . ,  m  —  1 

r  =  9  +  1,  • 

. . . ,  m 

*  =  1,2,... 

,«P 

T-H 

II 

.  ,nr. 

All  feasible  solutions  to  BLP  satisfy  BQP: 

Since  all  variables  are  binary,  any  feasible  solution  of  Equation  3.8  can 
be  adjusted  to  satisfy 


xipjq  —  xjqkri 


P  =  1,2,.. 

. ,  m  —  2 

9  =  P+  1,. 

. . ,  m  —  1 

r  =  9+  I,- 

. .  ,m 

*  =  1,2,... 

,np 

3  =  1,2,.. 

■  ,  nq 

Ar  =  1,2,.. 

•  ,  fir 

(3.19) 


such  that  there  is  no  change  in  the  objective  function  (Equation  3.2).  Therefore, 
from  Equation  3.10, 


"p 


'^/xtpjqxjqkr  ^  1, 

P=  1,2, ... 

. ,  m  —  2 

,=1 

9  =  p  +  1, . 

. . ,  m  — 

r  =  9  +  1,. 

. . ,  m 

j  =  1,2,... 

, ,  n, 

*  =  1,2,.. 

.  ,nr 

Up  Tip 

Tip 

^ '2,xipjqxjqkr  ~  5ZX»PJ92 

=  X^'PM’ 

P  =  1,2,.. 

. ,  m  —  2 

t=l  i=l 

i=l 

9  =  P+  I,- 

. . ,  m  — 

r  =  9+  1,. 

. . ,  m 

j  =  1,2,... 

.  ,n, 

*=  1,2,.. 

. ,  rip, 

(3.20) 


(3.21) 


showing  Equation  3.4  to  be  satisfied.  Applying  this  same  technique  to  Equations 
3.11  and  3.12  proves  Equations  3.5,  3.6,  and  3.15  to  be  satisfied.  And,  combining 
Equations  3.4  and  3.15  with  Equation  3.19  shows  Equation  3.3  to  be  satisfied,  as 
well.  Therefore,  the  two  problem  formulations  are  equivalent.  □ 
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3.2  Problem  Reduction 

One  approach  to  solving  the  track  initiation  problem  would  be  to  form 
all  possible  triples.  However,  the  computational  load  for  large  problems  would  be 
extensive,  regardless  of  the  method  used  to  solve  the  binary  linear  program.  To 
help  visualize  the  size  of  such  a  problem,  consider  a  track  initiation  buffer  of  m 
frames  with  the  number  of  unassigned  observations  in  frame  p  equal  to  np.  Then 
the  number  of  possible  triples  is 

m  —  2  m  —  1  m 

Y  Y  Ynpn<inr »  (3-22) 

p= 1  ?>P  r>9 

which  has  a  worst-case  complexity  of  0(m3n3),  where  n  is  the  number  of  actual 
targets.  For  a  100-target  scenario  with  a  buffer  of  seven  observation  frames, 
such  an  explicit  enumeration  would  involve  35,000,000  combinations.  In  practice, 
however,  most  of  the  unassigned  observations  will  be  assigned  to  clusters  within 
the  first  three  or  four  observation  frames.  This  conclusion  results  from  applying 
the  probability  of  a  successful  track  initiation, 

m  (m 

^(successful  initiation)  =  y](  . 

.=3  V  1 

where  p  is  the  probability  of  a  successful  detection  and  q  —  1  —  p,  to  various 
buffer  sizes  and  probabilities  of  missing  data.  Results  for  100  targets  are  shown 
in  Table  3.2.  But  even  under  these  conditions,  startup  could  require  in  excess 
of  1,000,000  combinations  be  examined.  Therefore,  the  track  initiation  process 


p  q 


(3.23) 


Buffer 

Size 

Probability  of  Missing 

0.05 

0.10 

0.20 

3 

85.74 

72.90 

51.20 

4 

98.60 

94.77 

81.92 

5 

99.88 

99.14 

94.21 

6 

99.99 

99.87 

98.30 

7 

100.00 

99.98 

99.53 

Table  3.2:  Expected  Number  of  Clusters  Initiated 


is  begun  by  first  reducing  the  number  of  allowable  arcs  based  upon  the  system 
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dynamics.  A  three-step  process  is  used  which  begins  by  forming  gates  around 
each  observation  to  limit  which  observations  can  be  linked  (Section  3.2.1).  For 
each  allowable  pair  of  linked  observations,  a  preliminary  estimate  of  the  specific 
energy  of  the  orbit  defined  by  that  pair  is  determined  and  evaluated  for  feasibility 
subject  to  the  system  dynamics  and  assumed  booster  characteristics  (Section 
3.2.2).  Then,  pairs  sharing  a  common  mid-point  are  linked  to  form  observation 
triples  and  initial  orbits  are  determined  (Section  3.3.1).  Each  triple  is  further 
evaluated  to  ensure  that  it  remains  feasible  in  terms  of  the  system  dynamics 
and  assumed  booster  characteristics  (energy).  Finally,  from  the  remaining  set  of 
feasible  triples,  the  binary  linear  program  defined  on  page  34  is  solved  (Section 
3.4). 


3.2.1  Single  Observation  Gating 

The  goal  of  this  section  is  to  determine  which  pairs  of  observations  may 
feasibly  be  linked  subject  to  system  dynamical  constraints.  To  accomplish  this 
goal,  estimates  of  pj,  pj,  aj,  and  £j  and  a  gate  for  each  estimate  are  calculated 
given  a  single  observation  at  time  t,  consisting  of  />,,  pi,  a,,  and  £,-.  These  estimates 
and  gates  are  formed  for  every  observation  such  that  tt  <  tQ  to  assess  links  with 
observations  in  each  observation  frame  where  t{  <  tj.  Actually,  this  process 
is  much  simpler  than  it  might  first  appear  since  the  targets  being  tracked  are 
assumed  to  be  acted  upon  only  by  the  gravity  of  a  spherical  earth.  Defining  the 
ECI  vector  to  the  target  as  r,  the  ECI  vector  to  the  sensor  as  r,,  and  the  range 
vector  from  the  sensor  to  the  target  as  p ,  then  the  three  vectors  are  related  by 
the  equation 

r  =  P  +  rs,  (3.24) 


as  seen  in  Figure  2.2. 

The  range  vector,  pt,  can  be  estimated  from  a  single  observation  as 


(cos  £,'  cos  Oti  ^ 
cos  £i  sin  a,  , 
sin£,  ) 


(3.25) 


where  p,  is  the  unit  range  vector.  Therefore,  the  sensor  position,  r,,  can  be 
estimated.  And,  since  the  only  acceleration  acting  on  either  the  target  or  the 
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sensor  is  assumed  to  be  due  to  the  earth’s  gravity, 


Ti  =  -—r,  and 

r 


r»  =  -rir> 
7 


(3.26) 


Now,  the  target’s  position  at  time  tj  can  be  estimated  through  the  use 
of  a  Taylor  series  expansion  as 

2 


Tj  =  r,  +  TiTji  r, 


+ 


(3.27) 


where  Tji  =  tj  —  <,.  While  not  enough  information  to  estimate  r*  is  available,  it 
is  possible  to  make  approximations  to  provide  reasonable  estimates  and  bounds 
for  those  estimates. 

Now,  by  definition 


where 


pj  =  II  pJ  =  v/pTpJ’ 


Pj  =  rJ  - 


(3.28) 


(3.29) 


Since  the  sensor’s  state  at  t}  is  known,  p3  can  be  expressed  (truncating  the  Taylor 
series  expansion  after  the  acceleration  term)  as 

1-2 


J • 


Pj  =  Fi  +  TiTji  +  r, 

=  4r,  +  TiTji  -  Tj„ 


where 


A=  1- 


Wi 
2  r,3 


Therefore, 

/>/  =  pj  ■  Pj  =  Mr.'  +  tiTji  -  Tj.)  ■  (/tr*  +  TiTji  -  Tj,), 
which,  when  expanded,  yields 

Pj2  =  A2( r,  •  r.)  +  2ArJl(rl  •  r.)  +  rJt-2(r,  •  r.)  -  2/t(r,  •  rja) 


-  2rJ,(r,  •  rJ3)  +  rJa 


=  Mr,)2  +  (v,TJt)2  -  2A(Ti  ■  r},)  +  rj,2  +  2Tji(ATi  -  Tj,)  ■  r, 


(3.30) 

(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3.35) 

(3.36) 


=  Mr,)2  +  ( Vi  Tj,  )2  -  2/t(r,  •  rJt)  +  rJS2  +  2tji\\At1  -  cos  /?.  (3.37) 


The  estimate  of  pj  is  completed  by  making  some  approximations  for  the 
terms  depending  on  r,.  Although  the  direction  of  the  velocity  vector  is  not  known, 
its  magnitude  can  be  estimated  in  the  problem  under  consideration.  This  result 
is  due  to  the  fact  that  the  target  is  assumed  to  be  operating  in  a  conservative 
force  field  and,  therefore,  its  energy  is  a  constant.  If  the  maximum  velocity  at 
burnout  is  assumed  to  be  Vb0  and  it  is  also  assumed  that  burnout  occurs  close 
to  the  earth’s  surface  where  r^o  is  approximately  the  earth’s  radius,  then  the 
maximum  specific  energy  is 

£  =  ~  —  (3.38) 

2 


and  the  velocity  at  any  subsequent  time,  such  as  is  given  by 


v'~f  V  + 


(3.39) 


The  only  term  still  not  known  is  the  value  of  /?,  the  angle  between 
(Ar,-  —  r,,)  and  r\.  The  angle  0  can  be  approximated,  however,  by  noting  that  in 
extreme  cases  0  =  0X  ±  02,  where  0X  is  the  angle  between  (Ar,-  -  ry,)  and  p{  and 
02  is  the  angle  between  p{  and  r,.  Therefore, 


Rx  =  cos-i  (**-*•)■  P* 
*  l!Ar,-r^|| 


(3.40) 


02  =  COS 


-t  r-  '  Pi 


(3.41) 


Although  0i  can  be  calculated  directly,  to  calculate  02  it  should  be 
noted  that  r,  •  p,  is  simply  the  target  velocity  along  the  line  of  sight,  u,|(,  and  that 


v,„  =  (Pi  +  r„)  •  Pi  =  pi  +  i-is  ■  Pi- 


Therefore, 


02  =  COS 


-1  P,  +  ■  Pr 


(3.42) 


(3.43) 


From  this  development,  the  remaining  uncertainty  in  the  calculation 
of  pj  results  from  not  knowing  the  relative  direction  of  the  components  of  the 
sensor  and  target  velocity  normal  to  the  line  of  sight.  The  extreme  values  of  0 
correspond  to  an  assumption  that  these  components  are  coplanar  and  in  either 


41 


the  same  or  opposite  directions.  Therefore,  the  larger  of  /?i  and  /?2  is  used  as  the 
nominal  value  and  the  smaller  to  determine  the  bounds  for  the  value  of  pj.  That 


P  =  max(/?i,  /?2)  ±  min^,  Pi)- 


Now,  to  estimate  pj, 


Pj  =  PjPj  =  l(Pj!!cosTi, 


(3-44) 


(3.45) 


where  7 y  is  the  angle  between  Pj  and  p-.  Since  pi  and  pj  are  known,  estimate 
can  be  completed  if  the  angle  between  p{  and  pJ  can  be  deter.  :  ed.  Designating 
this  angle  as  0,  the  distance  between  p,  and  pj  is 


tip  —  Jpi2  +  pj 2  -  2pipj  cos  0. 


(3.46) 


However,  the  angle  7 is  approximately  the  angle  between  6p  and  pj  and  since 


Pi2  =  fy2  +  Pj2  ~  2 tippj  cos  7 j. 


cos  7 j 


tip2  +  pj 2  -  pj2 
2  tippj 

Pj 2  -  PiPj  COS  0 


(3.47) 


(3.48) 


(3.49) 


Approximating  pj  as  ( tip/Tji ), 


Pj2  ~  PiPj  cos  0 

Pj  — 

PjTj> 


and  0  can  be  found  by  noting  that 


sin  0  =  — , 
Pj 


(3.50) 


(3.51) 


where  d±  is  the  distance  traveled  perpendicular  to  the  line  of  sight  between  tt 
and  tr  Therefore, 

„  vV  - 


COS0  — 


(3.52) 


ifl 

tt 


fflsa 
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Unfortunately,  the  azimuth  and  elevation  at  time  tj  cannot  be  estimated 
using  a  single  observation  since  the  target  velocity  vector  cannot  be  fully  esti¬ 
mated.  Therefore,  it  is  assumed  that  the  new  azimuth  and  elevation,  aj  and 
£j,  are  in  the  neighborhood  of  the  previous  azimuth  and  elevation,  cr,  and  e, 
(i.e.,  a.j  =  a,  and  £j  =  £,),  and  only  the  gates  for  these  measurements  will  be 
determined. 

The  maximum  angular  change,  6max,  will  have  components  in  both  az¬ 
imuth  and  elevation.  However,  since  there  is  no  way  of  knowing  exactly  which 
direction  the  target  is  moving  0max  is  used  as  the  gate  for  both  the  azimuth  and 
the  elevation,  resulting  in  a  somewhat  conservative  estimate.  Caution  must  also 
be  used  in  calculating  6a  when  £  «  ±|  since  a  small  change  of  9max  will  yield 
corresponding  large  changes  in  6a.  In  fact,  the  relationship  for  the  maximum 
change  in  a  for  a  given  0max  is 


and 


6a  =  (3.63) 

COS  £i 

6£  =  emax.  (3.64) 


Each  observation  in  an  observation  frame  is  assessed  in  this  manner  to 
determine  an  estimate  and  gate  for  each  observation  in  the  succeeding  observa¬ 
tion  frames.  Once  this  process  is  completed,  each  observation  in  the  succeeding 
observation  frames  is  compared  to  the  estimate  and  gate  for  that  frame.  This 
process  is  0(m2n2)).  If  the  observation  satisfies  the  gating  criteria,  the  pairing 
is  considered  for  addition  to  a  list  of  allowable  pairs. 


3.2.2  Dual  Observation  Gating 

Before  being  admitted  to  the  list  of  allowable  pairs,  a  dual  observation 
estimate  is  formed  to  determine  the  energy  of  the  orbit  defined  by  that  pair. 
Since  there  is  some  assumed  maximum  energy  which  can  be  achieved  based  upon 
the  known  (or  inferred)  booster  characteristics,  any  orbit  satisfying  this  energy 
restriction  is  considered  allowable  and  the  pair  is  added  to  the  list.  In  the  em¬ 
pirical  results  to  be  presented  later,  only  links  between  two  observations  of  the 
same  target  are  likely,  and  the  number  of  pairs  formed  appears  to  be  0(n). 


•  J 
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Because  only  the  energy  estimate  is  required  and  it  is  constant  for  a 
given  orbit,  the  estimates  for  rm  and  rm  can  be  determined  for  any  time  tm.  In 
this  section,  therefore,  two  complete  observations  of  range,  range  rate,  azimuth, 
and  elevation  at  times  ti  and  f2  will  be  used  to  determine  rm  and  rm  so  the  energy 
and  its  error  may  be  estimated. 

Recalling 


(cos  e,  cos  a, 
cos  e,  sin  a{ 
sin  £,■ 


i  =  1,2 


(3.65) 


r,  =  Pi  +  r.s, 


(3.66) 


The  range  vectors  iq  and  r2  can  be  expressed  using  Taylor  series  expansions  at 
some  common  epoch  tm  as 


Since 


1..  2 

**m  T  rtn T,m  T  2  ^w»7»m  T 


H 

*  = 


i  =  1,2. 


(3.67) 


(3.68) 


and  the  acceleration  vector  at  time  tm  can  be  expressed  using  Lagrange  interpo¬ 
lation  as 


..  Tm2 ..  Tm  1  .. 

rm  =  - r;  +  — r2 

r12  T21 

7*2  to  ..  7j  m  „ 

=  — rj - r2 

721  T21 

a  system  of  two  equations  with  two  unknowns  results 


m  \  I  4  m 


1  r2m  J  \  r, 


which,  when  solved  for  rm  and  rm  yields 

/  ~2m  Tin 


T21  T2i 

J_  J_ 

721  7"2i 


rmT2m 


(3.69) 


(3.70) 


(3.71) 


(3.72) 


'H,  \ JX  •«>  ’w-*  *.*>  J>k  L>  ->  W*>  V> 


From  this  estimate  of  the  state,  the  specific  energy  of  the  trajectory 
between  these  two  observations  can  be  estimated  according  to 


F  Vm  V 

m~  2  rm ' 


(3.73) 


To  determine  the  error  bounds  on  this  estimate  of  £m,  the  same  approach  as  used 
to  determine  S0  in  Section  2.1.4  is  applied.  In  this  case, 


Sm  =  Sm(  Sm) 


Sm  —  Sm(fl), 


(3.74) 


(3.75) 


where  ft  is  the  observation  set  required  for  the  state  estimate  and  may  be  written 


T 

=  (pi,  pi,  Oil  £«>  Sj,,  ti) 


(3.76) 


For  the  sake  of  analysis,  the  sensor  states,  Sis,  and  the  times  of  the 
observations,  ti,  are  again  assumed  to  be  known  exactly.  The  only  uncertainty 
is  associated  with  the  measurement  of  the  observation  attributes  pi ,  pi,  a,,  and 
£i.  It  is  further  assumed  that  each  attribute  type  is  independently  normally 
distributed  and  homoscedastic.  That  is 


Pi  ~  N(pp„ap), 
Pi  ~ 

oti  ~  N(p.a>,cra), 
and  £i  ~  N(pCi,<re). 


(3.77) 

(3.78) 

(3.79) 

(3.80) 


Then,  applying  a  Taylor  series  expansion  and  retaining  only  first-order 
terms,  it  can  be  shown  that 


where 


P  m  —  J  m  Pd  m  ^  , 


j  _  asj  _  dS2 

n=Hn  dil  n 


(3.81) 


(3.82) 
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(3.83) 


that 

where 


Knowing  Pm,  this  same  approach  can  be  applied  to  £m(Sm )  to  show 
<rf2  =  EmPmEmT,  (3.84) 


F  = 

m  dsm  sm=s* ’ 


For  this  specific  estimate,  where 

„  _  r2m  /  p  \  Tlm  (,  jj_  \ 

—  (1  o^lm^2m  )  1  1  )  ^2i 

r2 1  \  riJ  /  r2i  \  r2J  / 

=  ("l  nhra^2ml  **2  \  ^1) 

r21  V  r23  /  r2i  V  r,3  / 


(3.85) 

(3.86) 

(3.87) 
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(3.88) 


where 


drm  _  r2m  [7  n  \  dry  3/2  (  drA  ' 

dpi  ~  r21  [1  ri3TlmT2-J  dpi  +  ri5ri-T2-  [r*  '  Fl 

=  ^  [(i  -  ^rlmr2m^  pj  +  ^■rlmr2m(r1  •  pjfi 


pL  =  0 

dpi 


(3.89) 

(3.90) 

(3.91) 
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\  3 fx  (  drA  ' 

lmT2mJ  dot\  +  r!5TlmT2m  y1  '  daj  r\ 


3  Am^2m 
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(3.94) 

(3.95) 

(3.96) 

(3.97) 

(3.98) 

(3.99) 

(3.100) 
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(3.104) 


(3.105) 
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The  resulting  estimates  of  Sm  and  are  now  used  to  assess 
bility  of  the  observation  pair.  If 


(3.117) 
the  feasi- 


£m  <  £max  -f  3cr £■ , 


(3.118) 


then  the  observation  pair  is  added  to  the  list  of  allowable  pairs. 
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3.3  Track  Selection 


Once  the  list  of  allowable  pairs  is  formed,  the  list  of  allowable  triples 
is  developed.  To  link  two  pairs  of  observations  to  form  a  triple,  each  pair  must 
share  a  common  observation.  Therefore,  to  form  a  link,  the  second  observation 
of  the  first  pair  must  be  the  same  as  the  first  observation  of  the  second  pair.  The 
list  of  pairs  not  ending  in  the  current  observation  frame  is  evaluated  for  linkage 
with  pairs  in  succeeding  observation  frames  sharing  a  common  mid-point  and  for 
each  of  these  triples  an  initial  estimate  is  computed.  Once  again,  the  energy  of 
the  orbit  is  evaluated  as  an  additional  check  to  determine  if  the  triple  is  feasible. 
If  it  is,  the  calculated  state  and  state  covariance  are  stored  and  the  triple  is  added 
to  the  list  of  allowable  triples. 

If  only  0(n )  triples  are  considered,  then  it  is  possible  to  limit  the  com¬ 
plexity  of  this  phase  to  0(n)  by  judicious  application  of  indexing.  This  is  due 
to  the  fact  that  the  list  of  allowable  pairs  need  be  traversed  only  once  and  since 
it  is  known  that  the  starting  observation  of  the  second  pair  is  the  same  as  the 
ending  observation  of  the  first  pair.  Knowing  the  index  of  the  second  pair  in  the 
list  of  allowable  pairs  and  recalling  the  assumption  that  only  one  pair  is  likely  for 
a  given  observation  permits  each  triple  to  be  formed  directly. 


3.3.1  Initial  Estimate 


Given  that  an  observation  triple  is  formed,  estimates  of  the  target  posi¬ 
tion  and  velocity,  rm  and  rm,  at  some  time  tm  <  t0  are  now  calculated  using  three 
observations,  along  with  Sm  and  as  as  was  done  in  Section  3.2.2.  If  the  calculated 
specific  energy  satisfies  the  energy  restriction  of  Equation  3.118,  the  observation 
triple  will  be  added  to  a  list  of  allowable  triples  for  subsequent  consideration  as  a 
new  track  when  the  binary  linear  program  of  Section  3.1  is  solved  in  Section  3.4. 

Proceeding  with  Laplace’s  method,  the  position  vector,  rm,  may  be 

written 

=  Pm Pra  4”  fm)i  (3.119) 


and  the  velocity  vector,  rm,  is  found  by  taking  the  derivative  of  Equation  3.119, 


giving 


—  Pm  Pm  4"  PmPm  "b  r ms  ■ 


(3.120) 
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The  only  term  not  directly  available  from  the  sensor  measurement  set  is  pm, 
the  time  rate  of  change  of  the  unit  range  vector.  This  term  can  be  estimated, 
however,  by  using  Lagrange  interpolation  to  find  an  expression  for  the  unit  range 


vector 


*  a  7’ml^’m3  a  ^ml^m2  a 

Pm  =  ——Pi  +  —  —Pa  +  —  ~--p3» 

r12T13  t21t23  T32T31 


which,  when  differentiated  with  respect  to  tm  yields 

A  Tm 2  +  Tm 3  .  Tmi  +  rm3  A  Tm\  +  Tm 2  „ 

Pm  =  - Pi  +  — - P2  +  - P3- 

r12r13  r21r23  t32t31 


Moreover,  if  fm  =  f2, 


a  Tjt  A  Tji  —  T21  A  (  7"21  A 

P2  — - Pi  H - P2  d - P3? 

T2lT31  T21  Tji  7j,T31 


(3.121) 


(3.122) 


(3.123) 


everything  needed  to  calculate  the  initial  state  estimate  is  available.  It  is  not 
necessary  for  tm  to  equal  one  of  the  discrete  observation  times,  although,  if  it  did 
not,  it  would  also  be  necessary  to  apply  Lagrange  interpolation  to  find  estimates 
Pm)  Pm?  Pm  1  and  rma.  Therefore,  it  is  simpler  computationally  to  select 

the  central  observation  time  as  the  reference  ~poch. 

Proceeding  as  in  Section  3.2.2,  the  final  feasibility  check  is  performed 
on  the  observation  triple  using  the  energy  restriction  of  Equation  3.118.  That  is, 
£2  and  0£  are  computed  by  forming 


where 
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Specifically, 
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dp2  _  7~21  &P3 

dc*3  TjiT3\  das 

dp2 _  721  dp3 
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Then,  the  specific  energy  of  the  observation  triple  is 


and  its  variance  is 


C  _  vrn‘  p 

m  2  rm 


cr£2  =  EmPmEmT.  (3.144) 

If  the  calculated  energy,  £m,  satisfies  Equation  3.118,  the  observation  triple  is 
added  to  the  list  of  allowable  triples  and  the  final  stage  of  the  track  initiation 
process  is  performed. 

3.4  Cluster  Initiation 

With  the  set  of  allowable  groups  of  observations  to  initiate  clusters 
determined,  the  groups  which  provide  the  “best”  overall  assignment  must  now 
be  chosen.  To  begin  with,  since  a  rather  thorough  examination  of  each  triple 
has  been  performed,  it  can  be  concluded  that  each  triple  satisfies  the  physical 
constraints  of  the  system  dynamics  and  assumed  booster  characteristics,  and 
therefore,  the  final  solution  should  maximize  the  number  of  triples  selected.  And 
maximization  of  the  total  number  of  tracks  initiated  is  ensured,  because  the  total 
specific  energy  is  being  minimized  and  the  specific  energy  of  a  ballistic  orbit 
is  negative.  Selection  is  still  subject  to  the  overall  assignment  being  feasible, 
however.  That  is,  no  observation  can  belong  to  more  than  one  triple. 

Therefore,  if  the  initial  list  of  triples  is  feasible,  the  assignment  is  com¬ 
plete  and  a  new  cluster  is  initiated  for  each  triple.  If  the  overall  assignment  is 
infeasible,  however,  one  or  more  triples  must  be  removed  until  the  assignment 
is  feasible.  Since  this  process  may  result  in  more  than  one  feasible  solution,  a 
means  of  discriminating  among  these  to  select  the  “best”  one  is  also  required. 

Using  the  specific  energy  of  the  estimated  orbits  as  the  discriminant  as 
developed  in  the  binary  linear  program  of  Section  3.1,  the  set  of  clusters  with 
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the  minimum  total  specific  energy  is  defined  to  be  the  “best”  solution.  Such 
a  discriminant  is  possible  because  targets  on  ballistic  trajectories  follow  mini¬ 
mum  energy  paths  through  the  geopotential.  Any  misgrouping  of  observations  is 
assumed  to  result  in  a  higher  energy  orbit. 


To  actually  determine  the  optimal  cluster  assignments,  the  initial  set 
of  observation  triples  is  evaluated  for  feasibility.  If  it  is  not  feasible,  a  branch- 
and-bound  procedure  [28]  is  initiated  which  performs  a  depth-first-search  for  the 
optimal  solution.  One  cluster  is  removed  from  consideration  at  each  level  and 
a  branch  is  fathomed  when  it  becomes  feasible.  The  feasible  solution  with  the 
lowest  total  specific  energy  is  designated  as  the  optimal  solution  and  a  cluster  is 
initiated  for  each  active  triple.  The  worst-case  complexity  for  this  last  phase  is 
0(2")  if  the  entire  branch-and-bound  tree  must  be  searched. 

Through  the  multi-step  elimination  process  the  original  binary  linear 
program,  of  size  0(m3n3),  is  seen  empirically  to  be  reduced  to  0(m2n2)  by  de¬ 
termining  which  triples  are  dynamically  feasible.  The  final  process  of  selecting 
the  optimal  solution  is  expected  to  require  0(2")  time  in  the  worst  case  for  either 
approach. 

Once  the  track  initiation  process  is  complete,  all  new  clusters  are  added 
to  the  current  list  of  clusters  to  be  updated  through  the  track  maintenance  process 
and  their  observations  are  removed  from  the  pool  of  unassigned  observations. 
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Chapter  4 


Analysis  of  Results 


4.1  Simulation  Scenarios 

To  demonstrate  the  effectiveness  of  the  method  described  in  the  last 
two  chapters,  five  scenarios  are  developed.  These  scenarios  are  based  upon  the 
following  assumptions: 

•  Launch  points  are  randomly  generated  in  the  area  60°-70°  East  longitude, 
50°-60°  North  latitude,  an  area  covering  known  Soviet  ICBM  fields, 

•  100  targets, 

•  Impact  points  are  randomly  generated  in  the  area  80°-120°  West  longitude, 
30°-50°  North  latitude,  an  area  covering  most  of  the  continental  United 
States, 

•  Each  missile  has  identical  characteristics,  quantified  by  an  instantaneous 
Qbo  value  of  0.90  at  launch, 

•  Launch  occurs  randomly  over  the  first  30  seconds  of  the  scenario,  and 

•  A  spherical  rotating  earth  with  no  atmosphere, 

•  Simulation  begins  at  04:25:07  UTC  16  July  1986. 

These  scenarios  assume  that  an  ICBM  attack  is  being  observed  from  a  single 
ICBM  field  which  normally  contains  50  100  ICBMs  geographically  dispersed  to 
prevent  multiple  ICBMs  from  being  destroyed  in  a  single  strike.  The  value  of 
Qbo  was  based  upon  data  for  existing  Soviet  ICBMs  for  which  a  ballistic  missile 
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defense  is  being  contemplated1.  From  [12],  knowing  the  maximum  free-flight 
range  angle,  'F, 


Qbo  — 


2sin(fl/2) 

1  -F  sin('P/2) 


(4.1) 


Each  scenario  uses  the  same  launch  points  and  launch  times  and  pairs 
the  launch  points  with  impact  points  from  a  fixed  set  of  100  impact  points.  The 
scenarios  differ  in  which  impact  points  are  assigned  to  a  given  launch  point.  The 
trajectories  typical  of  these  scenarios  are  illustrated  in  Figure  4.1  which  shows  a 
polar  view  of  the  complete  trajectories  of  all  100  targets  used  in  Scenario  1. 

The  sensors  viewing  the  launch  scenarios  are  placed  in  modified  Molniya 
orbits  which  allow  all  targets  to  remain  in  view  of  the  sensors  for  the  duration 
of  each  of  the  five  scenarios.  The  specific  orbital  elements  used  are  given  in 
Table  4.1. 


Sensor 

Incli¬ 

nation 

Ascending 

Node 

Argument 
of  Perigee 

Mean 

Anomaly 

Eccen¬ 

tricity 

Mean 

Motion 

1 

85.0000 

70.0000 

281.3211 

132.4823 

0.3689374 

8.065369 

2 

85.0000 

140.0000 

281.3211 

132.4823 

0.3689374 

8.065369 

3 

85.0000 

210.0000 

281.3211 

132.4823 

0.3689374 

8.065369 

4 

85.0000-1 

350.0000 

281.3211 

132.4823 

0.3689374 

8.065369 

Table  4.1:  Sensor  Orbital  Elements 


Data  from  each  sensor  is  provided  at  five-second  intervals  for  the  first 
six  minutes  of  each  attack  scenario.  This  six-minute  period  covers  the  boost 
phase  and  post-boost  phase  prior  to  reentry  vehicle  deployment.  This  period  is 
of  particular  interest  in  an  SDI  scenario  because  it  is  easier  to  destroy  the  targets 
before  reentry  vehicle  deployment  due  to  the  smaller  number  and  larger  sizes  of 
the  targets. 

Figures  4.2  through  4.5  illustrate  the  complexity  of  the  problem.  These 
figures  give  a  good  indication  of  the  high  densities  of  the  targets  being  tracked. 
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1  Aviation  Week  &  Space  Technology ,  March  14,  1988,  page  153. 


They  also  show  large  numbers  of  target  crossings  and  similar  tracks  which  typi¬ 
cally  cause  the  most  difficulties  for  multi-target  tracking  algorithms.  Each  figure 
presents  the  view  from  one  of  the  four  sensors  for  the  first  six  minutes  of  the 
ICBM  attack  in  Scenario  1.  All  100  targets  are  present  in  each  view.  At  the 
bottom  of  each  figure,  the  information  pertaining  to  the  sensor  vantage  point  is 
displayed,  showing  the  latitude  and  longitude  of  the  satellite  sub-point  along  with 
the  altitude  of  the  satellite  above  the  earth’s  surface  and  the  simulation  time  of 
that  position. 


76.0  N  90.6  E  728/. 5  km  0.00  sec 

Figure  4.2:  View  from  Sensor  1 — Scenario  1 


76.0  N  129.4  W  7287.5  km  0.00  sec 

Figure  4.4:  View  from  Sensor  3 — Scenario  1 


Four  variance  sets  and  four  levels  of  missing  data  were  considered  during 
the  course  of  this  investigation  to  examine  the  effects  of  high  and  low  quality  data. 
The  scenario  variances  are  listed  in  Table  4.2.  The  four  levels  of  missing  data  are 
0,  5,  10,  and  20  percent. 


Variance 

Set 

Variances 

Range 

(meters)2 

Range  Rate 
(meters/second)2 

Azimuth 

(radians)2 

Elevation 

(radians)2 

1 

10.0 

1.0 

io-7 

IO"7 

2 

10.0 

1.0 

10-8 

00 

1 

o 

3 

10.0 

1.0 

10"9 

10"9 

4 

10.0 

1.0 

o 

rH 

1 

o 

t— H 

IO"10 

Table  4.2:  Scenario  Variances 

To  evaluate  the  performance  of  the  temporal  clustering  algorithm,  two 
sets  of  variances  and  levels  of  missing  data  are  used.  The  first  case,  representing 
expected  values  for  the  measurement  variances  and  levels  of  missing  data,  uses 
Variance  Set  3  and  5  percent  missing  data,  while  the  second  case ,  representing 
worst-case  values,  uses  Variance  Set  1  and  20  percent  missing  data. 

4.2  Discussion  of  Results 

The  temporal  clustering  procedure  described  in  the  previous  two  chap¬ 
ters  is  implemented  in  Pascal  on  the  Cray  X-MP/24  at  The  University  of  Texas 
at  Austin’s  Center  for  High  Performance  Computing.  A  listing  of  the  code  used 
is  included  in  Appendix  B.  Before  the  code  can  be  run,  however,  determinations 
of  the  size  of  the  track  initiation  buffer  and  the  maximum  number  of  missing  ob¬ 
servations  for  track  termination  are  necessary.  Based  upon  the  results  presented 
in  Table  3.1  on  page  32  and  Table  2.1  on  page  27,  a  track  initiation  buffer  of  seven 
observation  frames  and  a  maximum  of  five  consecutive  missing  observations  was 
selected  as  a  termination  criterion  for  the  maximum  value  of  20  percent  missing 
data,  since  the  expected  number  of  track  initiation  failures  and  false  terminations 
over  the  life  of  each  of  the  scenarios  investigated  is  less  than  one. 


Tables  4.3  and  4.4  are  summaries  of  the  results  of  the  twenty  simulation 
runs  (five  scenarios  with  four  sensors  each)  used  to  investigate  the  performance 
of  the  temporal  clustering  algorithm  for  both  the  expected  and  worst-case  values 
of  measurement  noise  and  missing  data. 


Scenario/ 

Sensor 

Perfect 

Unassigned 

Observations 

Errors 

Observations 

Clusters 

Targets 

Termination 

Gating 

1/1 

100 

100 

0 

0 

0 

6582/6944 

1/2 

100 

99 

1 

0 

0 

6588/6944 

1/3 

100 

100 

0 

0 

0 

6591/6944 

1/4 

100 

100 

0 

0 

0 

6647/6944 

2/1 

100 

99 

1 

0 

0 

6582/6944 

2/2 

100 

99 

1 

0 

0 

6588/6944 

2/3 

100 

100 

0 

0 

0 

6591/6944 

2/4 

100 

100 

0 

0 

0 

6647/6944 

3/1 

100 

99 

1 

0 

0 

6582/6944 

3/2 

100 

99 

1 

0 

0 

6588/6944 

3/3 

100 

100 

0 

0 

0 

6591/6944 

3/4 

100 

100 

0 

0 

0 

6647/6944 

4/1 

100 

99 

1 

0 

0 

6582/6944 

4/2 

100 

99 

1 

0 

0 

6588/6944 

4/3 

100 

100 

0 

0 

1B9HI 

6591/6944 

4/4 

100 

100 

HHKIHi 

0 

6647/6944 

5/1 

100 

0 

0 

6582/6944 

5/2 

100 

99 

i 

0 

0 

6588/6944 

5/3 

100 

100 

0 

0 

0 

6591/6944 

5/4 

100 

100 

0 

0 

0 

6647/6944 

2000 

1992 

8 

0 

0 

132040/138880 

Table  4.3:  Summary  of  Results — Nominal  Values 


The  first  column  of  the  table  indicates  the  scenario  number  and  sensor 
number  of  the  observations.  The  next  two  columns  show  the  number  of  perfect 
clusters  and  perfectly  clustered  targets  for  each  run.  A  cluster  is  considered 
perfect  if  all  of  its  observations  are  from  a  single  target.  No  imperfect  clusters 
were  encountered  in  any  of  the  runs  examined.  A  target  is  considered  perfectly 
clustered  if  all  of  its  observations  are  contained  in  a  single  cluster. 

On  average,  each  run  in  the  nominal  case  contained  less  than  one  imper¬ 
fectly  clustered  target,  resulting  from  a  single  observation  being  excluded  from 


LHTVX  VXATXVX  VKMK  IT*  1 


Scenario/ 

Sensor 


Perfect 

Clusters  |  Targets" 


Unassigned  _ Errors _ 

Observations  Termination  Gating  Observations 

1  1  1111  5556/6944 

0  0  1  5568/6944 

_ _ 0 _ 1 _ 0  5536/6944 

2  0  1  5561/6944 

__  1 _ 1 _ 1  5556/6944 

0  02  5568/6944 

0  1  0  5536/6944 

2  0  1  5561/6944 

1  1  1  5556/6944 

0  0  1  5568/6944 

0 _ 1 _ 0  5536/6944 

2  _ 0 _ 1  5561/6944 

1  1  1  5556/6944 

0 _ 0 _ 2  5568/6944 

0 _ 1 _ 0  5536/6944 

2  _ 0  1  5561/6944 

2 _ 1  1  5556/6944 

0 _ 0 _ 1  5568/6944 

0  1  0  5536/6944 

2  |  0  1  1  |  5561/6944  ~ 

16  10  |  17  |  111105/138880 


Table  4.4:  Summary  of  Results — Extreme  Values 

a  gate  either  during  the  track  maintenance  or  track  initiation  phase.  The  fre¬ 
quency  of  occurrence  is  shown  in  column  four  of  Table  4.3.  However,  all  clusters 
developed  were  perfect.  Considering  the  large  number  of  target  measurements 
observed  (132,040)  in  the  twenty  cases  evaluated,  losing  only  eight  observations 
is  an  exceptional  performance. 

Using  the  worst-case  values,  each  run  averages  one  split  cluster  and 
two  improperly  clustered  targets,  still  quite  good  performance.  These  failures 
result  from  one  of  three  causes  and  their  frequencies  are  shown  in  columns  4-6  of 
Table  4.4.  The  first  cause  is  the  result  of  an  observation  not  being  assigned  to  any 
cluster.  In  five  cases,  two  observations  were  lost  during  the  track  initiation  process 
because  of  a  failure  to  obtain  three  observations  in  the  seven-frame  track  initiation 
buffer.  Each  of  these  failures  could  have  been  avoided  by  simply  expanding 
the  size  of  the  buffer  but  the  increased  computational  complexity  and  storage 
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overhead  may  not  merit  such  a  change.  This  result  is  not  unexpected  based  on 
the  numbers  in  Table  3.1. 


Additionally,  one  observation  of  a  target  was  not  included  in  any  cluster 
because  the  problem  reduction  phase  of  the  track  initiation  process  incorrectly 
eliminated  it  from  consideration  in  a  feasible  track  based  upon  system  dynamics. 
The  remaining  five  unassigned  observations  were  the  result  of  correct  assignments 
being  discarded  during  the  track  maintenance  process  because  they  failed  the 
gating  criteria. 

The  other  reasons  for  incorrectly  clustered  targets  are  due  to  track  ter¬ 
mination  as  the  result  of  a  termination  or  gating  error.  A  termination  error  occurs 
when  a  track  is  terminated  because  more  than  the  maximum  number  of  missing 
observations  occurred.  Again,  this  type  of  failure  could  be  avoided  by  increasing 
the  maximum  number  of  missing  observations  with  a  corresponding  trade-off  in 
computational  cost.  Increasing  the  buffer  size  to  eight  observation  frames  would, 
therefore,  have  eliminated  ten  tracking  errors. 

A  gating  error  occurs  when  a  track  is  terminated  because  the  estimate 
of  a  target’s  state  is  sufficiently  far  from  the  true  state  to  prevent  the  correct 
observations  from  falling  within  the  observation  gates.  All  seventeen  gating  errors 
occurred  immediately  after  track  initiation  during  the  EKF  stabilization  phase. 
These  gating  errors  could  be  reduced  (or  effectively  eliminated)  through  either 
reduced  measurement  noise  (improved  sensor  characteristics)  or  improved  orbit 
determination  procedures.  This  conclusion  is  supported  by  the  total  lack  of  gating 
errors  in  the  cases  using  nominal  values  of  measurement  noise  and  missing  data. 

And,  as  expected,  lower  measurement  noise  results  in  better  estimates  of 
the  target  state.  The  effect  of  various  levels  of  measurement  noise  on  the  estimates 
of  target  position  and  velocity  are  illustrated  in  Figures  4.6  and  4.7.  Each  figure 
shows  the  difference  between  the  true  and  estimated  position  or  velocity  for  the 
four  levels  of  measurement  noise  considered.  All  data  is  for  Scenario  1,  Sensor  1, 
Cluster  1. 

As  a  final  test  of  the  robustness  of  the  temporal  clustering  algorithms, 
biases  were  introduced  into  the  sensor  attitude  and,  independently,  in  the  sensor 
clock  and  applied  to  the  worst-case  runs.  An  attitude  bias  is  possible  due  to 
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Figure  4.6:  Position  Error  Due  to  Measurement  Noise 
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Figure  4.7:  Velocity  Error  Due  to  Measurement  Noise 


precessional  and  nutational  effects  of  the  sensor  platform  during  the  course  of 
its  orbit.  A  constant  bias  of  up  to  60  arc  seconds  results  in  only  two  additional 
unassigned  observations  in  the  twenty  runs  considered,  both  failures  in  the  track 
initiation  process.  These  biases  are  much  larger  than  one  would  expect  in  prac¬ 
tice,  since  the  slowly  varying  effects  would  be  expected  to  be  small  and  daily 
observation  of  the  sensor  platform  attitude  by  ground  personnel  would  allow  the 
development  of  models  to  remove  most  of  the  remaining  bias. 

A  clock  bias  would  result  in  the  sensor  thinking  it  was  in  a  position  dif¬ 
ferent  from  that  indicated  by  its  onboard  ephemeris.  Typically,  with  the  atomic 
clocks  available  on  most  satellites  these  biases  would  be  on  the  order  of  mi¬ 
croseconds  or,  at  worst,  milliseconds.  Yet,  the  temporal  clustering  algorithms 
performed  exactly  as  shown  in  Table  4.4  with  biases  of  up  to  0.5  seconds,  far 
above  any  bias  which  could  reasonably  be  expected. 

While  the  temporal  clustering  algorithm  performs  admirably  even  un¬ 
der  these  adverse  conditions  of  attitude  and  clock  biases,  the  effects  on  the  state 
estimate  and  state  covariance  are  considerably  degraded,  as  can  be  seen  in  Fig¬ 
ures  4.8  through  4.11.  Each  figure  shows  the  difference  between  the  true  and 
estimated  position  or  velocity  for  various  levels  of  bias.  All  data  is  for  Scenario  1, 
Sensor  1,  Cluster  1.  Obviously,  while  the  algorithms  perform  well  under  these  cir¬ 
cumstances,  it  is  highly  desirable  to  reduce  the  sensor  biases  as  much  as  possible 
to  attain  the  high  degree  of  accuracy  necessary  to  direct  a  proper  response. 
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Figure  4.9:  Velocity  Error  Due  to  Attitude  Bias 


Chapter  5 


Conclusions 


The  temporal  clustering  algorithm  developed  in  Chapters  2  and  3  does 
indeed  accomplish  the  initial  objectives  that  it 

•  Perform  both  track  initiation  and  track  maintenance  and 

•  Permit  processing  of  data  in  real  time  while  minimizing 

—  Computational  complexity  and 
-  Data  storage  requirements. 

Through  prudent  application  of  existing  solution  techniques  in  astrodyn'-.nics, 
mathematical  programming,  numerical  analysis,  and  statistical  estimation,  an 
integrated  solution  is  developed  which  is  not  only  capable  of  performing  both 
track  initiation  and  track  maintenance,  but  also  improves  on  previous  work  in 
the  fields  of  multi-target  tracking  and  clustering  to  effectively  track  large  numbers 
of  targets  in  real  time. 

In  fact,  the  duty  cycle  for  each  run  (the  ratio  of  the  total  execution 
time  to  the  elapsed  scenario  time)  is  only  4.4  percent.  Execution  time  for  all  runs 
averages  15.9  CPL  seconds  per  run.  The  portion  of  that  time  spent  in  each  major 
procedure  discussed  in  Figure  1.2  is  provided  in  Table  5.1.  Considering  that  only 
17.5  percent  of  of  the  total  execution  time  is  taken  up  by  algorithms  of  complexity 
worse  than  0(n)  and  76.7  percent  is  devoted  to  algorithms  which  are  vectorizable 
and  capable  of  being  run  in  parallel,  the  timing  results  are  quite  impressive. 
In  fact,  empirical  results  from  runs  with  20,  50,  and  100  targets  showed  the 
complexity  of  Perform  Cluster  Assignments  and  Perform  7'rack  Initiation  to  be 
only  0[n2)  and  0(n ),  respectively. 


Procedure 

Percent  Total 
Execution  Time 

Complexity 

1)  Read  Observation  Frames 

5.8 

0(k) 

2)  Forecast  Existing  Clusters 

47.5 

0(n) 

3)  Calculate  Assignment  Costs 

6.5 

0(n2) 

4)  Perform  Cluster  Assignments 

8.5 

0(n3) 

5)  Update  Clusters 

29.2 

0(n) 

6)  Perform  Track  Initiation 

2.5 

0(2") 

Table  5.1:  Timing  Results 

It  should  be  cautioned,  however,  that  these  timing  results  apply  only 
to  the  specific  case  investigated  in  this  study.  Any  combination  of  measurement 
types  and  system  dynamics  which  does  not  permit  a  significant  reduction  in  the 
size  of  the  binary  linear  program  (as  was  done  in  Section  3.2)  may  not  be  able 
to  satisfy  the  requirement  that  the  problem  be  solved  in  real  time.  In  such 
circumstances,  more  sophisticated  methods  of  solving  the  binary  linear  program 
may  be  implemented  to  reduce  the  computational  complexity,  although  there  is 
no  guarantee  that  these  improved  techniques  will  allow  real  time  processing  in 
all  cases. 

The  limitations  due  to  computational  complexity  and  data  storage  re¬ 
quirements  are  reduced  through  the  judicious  application  of  existing  algorithms 
for  filtering  (the  Extended  Kalman  filter),  assignment  (the  Hungarian  method), 
and  quadratic  programming  (branch-and-bound)  while  taking  full  advantage  of 
the  temporal  component  of  the  data  and  system  dynamics.  And,  as  seen  in  Sec¬ 
tion  2.3  and  at  the  beginning  of  Chapter  3,  additional  reductions  in  complexity 
and  data  storage  requirements  are  possible  if  the  missing  data  rate  is  kept  small. 

A  most  remarkable  feature  of  the  temporal  clustering  algorithm  is  its 
ability  to  function  well  when  faced  with  low  data  rates  and  high  levels  of  both 
missing  data  and  measurement  noise.  Even  in  the  runs  examined  using  the  worst- 
case  variances  (Variance  Set  1)  and  20  percent  missing  data,  the  temporal  clus¬ 
tering  algorithm  successfully  clustered  nearly  100  percent  of  the  over  110,000 
observations  available.  And  for  tracks  with  four  or  more  observations  all  but  10 


out  of  2000  targets  were  tracked  correctly  throughout  the  scenarios,  and  those 
were  the  result  of  exceeding  the  size  of  the  track  termination  buffer — something 
which  can  be  avoided  by  simply  increasing  the  size  of  the  buffer  by  one  frame. 
And  even  the  addition  of  biases  in  the  sensor  attitude  or  clock  to  levels  well  in 
excess  of  those  which  can  be  reasonably  expected  failed  to  significantly  affect  the 
performance  of  the  algorithm. 

5.1  Future  Research 

There  remain  areas  in  which  the  temporal  clustering  algorithm  can  be 
improved  or  the  scope  of  its  application  broadened.  More  attention  can  be  applied 
to  improving  the  overall  efficiency  of  the  tracking  process  through  the  applica¬ 
tion  of  state-of-the-art  filters  such  as  those  discussed  by  Kaminski,  Bryson,  and 
Schmidt  [35]  and  Verhaegen  and  Van  Dooren  [57]  or  simplification  of  the  filters 
through  the  application  of  constant  gain  Kalman  filters  as  discussed  by  Blackman 

[13]. 

In  addition,  significant  advantage  can  be  gained  by  exploiting  the  par¬ 
allel  structure  of  many  existing  computers  and  the  application  of  pipelining  as 
discussed  by  Allen,  Kurien,  and  Washburn  [1].  For  example,  considerable  im¬ 
provement  in  processing  time  can  be  achieved  by  developing  a  parallel  structure 
capable  of  independently  tracking  each  target,  especially  since  massively  parallel 
architectures  with  65,536  processors  exist  today.  This  is  particularly  true  consid¬ 
ering  that  over  75  percent  of  the  total  execution  time  for  the  temporal  clustering 
procedure  is  used  by  the  forecasting  and  state  update  algorithms.  These  func¬ 
tions  can  easily  be  performed  on  separate  parallel  processors  for  each  cluster. 
Use  of  vector  processing  is  also  helpful  when  integrating  the  large  state  and  state 
covariance  vectors  in  the  forecasting  algorithm. 

And  while  the  case  presented  assumes  a  spherical  earth  with  no  drag 
and  no  thrust,  the  method  can  be  readily  extended  to  cases  using  higher  order 
gravitational  potentials  and  atmospheric  drag  by  simply  reformulating  the  fil¬ 
ter  and  track  initiation  gating  process  to  specifically  account  for  these  effects. 
And,  depending  upon  the  specific  system  dynamics,  it  may  also  be  necessary  to 
choose  another  cost  coefficient  for  the  objective  function  of  the  binary  linear  pro- 


gram.  Reformulation  to  account  for  other  measurement  combinations,  such  as 
the  angles-only  case,  is  possible,  although  preliminary  investigations  have  shown 
this  case  to  be  somewhat  more  complex  due  to  marginal  observability.  Investiga¬ 
tions  of  vehicle  thrust  and  maneuvering  are  certainly  possible,  as  well,  although 
the  reformulation  will  likely  be  considerably  more  difficult  and  require  the  use  of 
adaptive  filtering  techniques.  Additional  investigations  of  algorithm  performance 
in  the  face  of  changing  measurement  covariance  due  to  sensor  degradation  or 
failure  are  also  possible  through  the  application  of  adaptive  filtering. 

And,  finally,  the  extension  of  this  temporal  clustering  approach  to  the 
broader  issue  of  multi-sensor  correlation  should  be  straightforward.  Although 
issues  of  distributed  processing  [21,22,23]  need  to  be  examined  in  detail,  it  appears 
that  once  a  cluster  is  established,  its  state  estimate  and  state  covariance  matrix 
can  be  transmitted  by  each  sensor  to  a  central  processor  for  correlation  with  data 
from  other  sensors.  Transmission  of  this  minimal  amount  of  data  significantly 
reduces  the  bandwidth  required  for  data  exchange  and  the  correlation  process 
can  then  apply  gating  and  assignment  procedures  quite  similar  to  those  used  in 
the  single- sensor  case. 

While  the  temporal  clustering  process  as  developed  here  is  specifically 
tailored  for  a  single  application,  prudent  modification  of  the  application-specific 
portions  should  allow  it  to  be  applied  to  other  ballistic  tracking  problems  or  even 
those  in  the  areas  and  fluid  dynamics  or  particle  physics. 
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Appendix  A 
Simulation  Design 


To  permit  realistic  assessment  of  the  approach  presented  here,  some 
means  of  providing  a  set  of  target  measurements  was  necessary.  To  accomplish 
this  objective,  a  data  generating  program  (GENDAT)  was  developed  by  the  au¬ 
thor  and  Stuart  H.  Smith.  A  flowchart  of  GENDAT  is  presented  in  Figure  A.l. 

GENDAT  was  designed  to  generate  both  the  target  and  sensor  states 
using  a  Runge-Kutta  4(5)  integrator  and  a  user-provided  force  model.  Input  to 
GENDAT  begins  by  determining  the  booster  characteristic  thrust,  Qb0 ,  where 


Qbo  ~ 


Vbo  r bo 


(A.l) 


and  whether  high  or  low  trajectories  are  to  be  used  for  the  targets.  Then,  the 
launch  time  and  the  launch  and  impact  points  for  each  target  are  input.  Finally, 
the  Keplerian  orbital  elements  for  each  sensor  are  input. 

An  initial  state  estimate  is  then  determined  for  each  target  at  its  spec¬ 
ified  launch  time  by  computing  the  trajectory  necessary  to  reach  the  assigned 
impact  point  from  the  designated  launch  point.  Initial  state  estimates  for  the 
sensors  are  also  determined  by  converting  the  orbital  elements  at  the  initial  sim¬ 
ulation  time. 

Once  all  the  initial  states  are  calculated,  GENDAT  begins  forming  ob¬ 
servations  of  each  target  from  each  sensor  while  the  target  is  in  view.  Currently, 
the  target  is  in  view  if  it  has  launched  and  not  yet  reached  its  impact  point.  No 
additional  considerations  such  as  sensor  field-of-view,  sensor  range,  or  obscura¬ 
tion  by  the  earth’s  limb  are  yet  implemented.  Observations  are  calculated  using 
the  transformation  described  in  Section  2.1.4. 
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Figure  A.l:  GENDAT  Flowchart 

For  each  measurement  attribute,  a  Gaussian  random  variable  is  com¬ 
puted  to  allow  for  expected  measurement  noise.  These  variables  are  N( 0, 1).  In 
addition,  a  single  uniform  random  variable  is  also  computed  for  use  in  simulating 
the  stochastic  nature  of  detecting  an  observation. 

The  target  and  sensor  states  are  output  to  separate  files,  the  former 
to  be  available  for  comparison  with  estimated  target  states,  the  latter  to  act 
as  the  satellite  ephemeris.  The  target  measurements  and  associated  random 
variables  are  output  to  a  third  file.  The  separation  of  the  true  measurements  from 
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the  measurement  noise  allows  the  temporal  clustering  algorithm  to  incorporate 
differing  measurement  variances  at  runtime. 

After  all  observations  are  formed  and  the  resulting  data  is  output,  the 
target  and  sensor  states  are  integrated  to  the  next  observation  time  and  the 
process  is  repeated  until  the  simulation  end  time  or  all  targets  have  reached  their 
destinations. 


Appendix  B 


Program  Listing 


Program  TClust er (tmdata , s sdata , t sdat a, t cldata , trudata , csdata , input , output ) ; 

{  Author:  TS  Kelso 

{  Original  Version:  26  January  1987  > 

{  Current  Revision:  26  Hay  1988  > 

{  Program  Description:  Program  performs  temporal  clustering  on  time 

successive  data  frames  using  an  assignment  from  last 
member  of  existing  clusters  to  observations.  Track 
initiation  performed  using  quadratic  program.  Assumes 
range,  range  rate,  azimuth,  and  elevation 
measurements.  ) 

(*#A+:R-  *) 


const 


clusters 

= 

201; 

{Maximum  allowable  clusters  +  1} 

max_pairs 

= 

400;  {Maximum  number  of  pairs  in  Perform.Cluster .Initiation) 

nest 

- 

6; 

{Elements  in  state  vector) 

block 

= 

3; 

{Axes  in  ECI  coordinate  system) 

nstack 

= 

42; 

{Elements  in  stacked  state  vector  =  nest  +  nest*2) 

nterms 

= 

8; 

{Maximum  number  of  terms  in  estimates) 

attributes 

= 

4; 

{Maximum  number  of  possible  measurement  attributes) 

bad 

= 

5.0; 

{Metric  value  for  bad  assignment  =  attributes  +  1) 

frame 1 

= 

-6; 

{Data  required  for  initial  estimate  =  framel.. 0) 

max_prop 

= 

5; 

{Maximum  number  of  propagation  intervals) 

max .missed 

= 

i; 

{Maximum  number  of  missed  gates  allowed) 

zero 

= 

1.0E-14; 

{Machine  epsilon  for  real  =  double) 

big 

= 

1.0E+14; 

mu 

= 

3.9860064E+14;  {Geocentric  gravitational  parameter,  m“3/s"2> 

small 

= 

1.0E-12; 

{RK78  integration  control  factor) 

Pi 

= 

3.1415926535897932; 

max _ energy 

= 

-3.437E7; 

{Maximum  specific  energy,  meters~2/second"2) 

type 

span 

=  framel. 

•  0; 

atr. vector  =  array  [0 .. attributes]  of  real; 

obs_vector  =  array  [1 .. attributes]  of  real; 

state_vector  =  array  [l..nest]  of  real; 
stacked.vector  =  array  [l..nstack]  of  real; 


=  airray  [1..2]  of  obs.vector; 

=  aurray  [1 .. clusters]  oi  boolean; 

=  array  [0. .clusters]  of  integer; 

=  array  [1  .clusters]  of  atr_vector; 

=  array  [1. . nest.l. .nterms]  of  real; 

=  array  [1. .nterms, 1. .nterms]  of  real; 

=  array  [1 . .nest, 1 . .nest]  of  read; 

=  array  [1 .. attributes , 1 . .nest]  of  real; 

=  array  [1. .nest, 1. .attributes]  of  real; 

=  array  [1. .attributes, 1. .attributes]  of  real; 
=  array  [spam]  of  state_vector; 

=  array  [spam]  of  obs_vector; 

=  array  [1. .clusters, 1. .clusters]  of  real; 


limits  =  array  [1..2] 

bvector  -  array  [1 . . clu 

ivector  =  airray  [0..clu 

frame  =  aurray  [1  .clu 

J_matrix  =  airray  [l..nes 

ER_matrix  =  array  [l..nte 

P_matrix  =  array  [l..nes 

H_matrix  =  array  [l..att 

K_matrix  =  array  [l..nes 

R_matrix  =  array  [l..att 

S_matrix  =  aurray  [spam] 

0_matrix  =  array  [spam] 

M_matrix  =  array  [1 . . clu 

states  =  record 

number  :  integer; 

time  :  real; 

values  :  state_vector ; 

end ;  {record} 
measures  =  record 

target, sensor  :  integer; 
time  :  real; 

obs.error  :  obs_vector; 
missing  :  real; 

end;  {record} 


sensor , sen.nr .max.target , 

next .target , last .target .time 

max.t ime , min.t ime , tbias , 

step, missing.f lag, missing.limit 

nr.unassigned.nr.obs , 

nr.clusters , nr .act i ve , nr.inact i ve 

frame.time 

nr .miss ing , convert , 

col.baais .row.basis 

status,  observation,  tairgets 

assigned 

R,bias 

spam.limits 

next 


Ident ity ,Q_k 
R.k 

metric,  aunetric 

Rvar 

attr 

est .gate 


boolean; 


integer; 


{End  of  input} 


array  [span]  of  integer; 
array  [span]  of  real; 

ivector; 

array  [spam]  of  ivector; 
array  [framel..i]  of  bvector; 
obs.vector ; 
limits ; 
atr.vector ; 

array  [1..59]  of  real; 

array  [spam]  of  state.vector ; 

P.matrix; 

R.matrix; 

M.matrix; 

array  [i . . 3]  of  ER.matrix; 

array  [spam]  of  frame; 

array  [1 . .clusters]  of  obs.vector; 


<  a 


Sx 

P 

Stan 

atr.limits 
current _obs 

current.sensor, current .target 
tmdata 

ssdata.tsdata 
tcldata, trudata, csdata 
times 

last , total 


array  [1 .. clusters]  of  state.vector; 
array  [1 .. clusters]  of  P .matrix; 
array  [span, 1. .clusters]  of  obs.vector; 
array  [framel..l]  of  limits; 
measures; 
states ; 

file  of  measures; 
file  of  states; 
text ; 

array  [0..9]  of  integer; 
array  [0..9]  of  real; 


{***  Timing  Functions  a*****************************************************} 

Function  Second  :  real;  FORTRAH; 

Procedure  Init .Times; 
var 

i  :  integer ; 
begin 

for  i  :=  1  to  9  do 
begin 

last[i]  :=  0.0; 
total [i]  :=  0.0; 
times [i]  :=  0; 
end;  {for  i> 

end;  {Procedure  Init .Times} 

Procedure  Start_Timer(arg  :  integer); 
begin 

last[arg]  :=  Second; 

end;  {Procedure  Start .Timer} 

Procedure  Stop_Timer(arg  :  integer); 
var 

elapsed  :  real; 
begin 

elapsed  :=  Second  -  last[arg]; 
total [arg]  :=  total [arg]  +  elapsed; 
times [arg]  :=  times [arg]  +  1; 
end;  {Procedure  Stop .Timer} 

Procedure  Report .Times (arg  :  integer); 
var 

i  :  integer; 
begin 

for  i  :=  1  to  arg  do 
begin 

Vrite(i : 2, ’ )  ■); 
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case  i  of 

1  :  WritelnC ’Input  Data’); 

2  :  WritelnC ’Forecast’ )  ; 

3  :  WritelnC ’Calculate  Metrics’); 

4  :  WritelnC ’Perform  Cluster  Assignment’); 

5  :  WritelnC ’Perform  Cluster  Initiation’); 

6  :  WritelnC ’Update  Estimates’); 
end;  {case} 

WritelnC’  Elapsed  time  =  ’ ,total[i] :7:4, 

’,  Average  time  =  * .total [i] /times [i] : 7:4, 

’,  Percentage  =  ’ , 100*total [i] /total [0]  : 4:1,  ’*/,’) ; 
end;  {for  i} 

WritelnC’  Total  time  =  ’ ,total[0] :7:4) ; 

Writeln; 

end;  {Procedure  Report .Times} 

{***  Global  routines  ******************************************< 

Function  IMinCargl ,arg2  :  integer)  :  integer; 
begin 

if  argl  <  arg2  then 
IMin  :=  argl 
else 

IMin  :=  arg2; 
end;  {Function  IMin} 

Function  IMaxCargl ,arg2  :  integer)  :  integer; 
begin 

if  argl  >  arg2  then 
IMax  :=  argl 
else 

I Max  : =  arg2 ; 
end;  {Function  IMax} 

Function  RMinCargl ,arg2  :  real)  :  real; 
begin 

if  argl  <  arg2  then 
RMin  :=  argl 
else 

RMin  :=  arg2; 
end;  {Function  RMin} 

Function  RMaxCargl ,arg2  :  real)  :  real; 
begin 

if  argl  >  arg2  then 
RMax  : =  argl 
else 

RMax  : =  arg2 ; 
end;  {Function  RMax} 
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Procedure  Increment (var  arg  :  integer); 
begin 

arg  :=  arg  +  1; 

end;  {Procedure  Increment} 

Procedure  Decrement (var  arg  :  integer); 
begin 

arg  :=  arg  -  1; 

end;  {Procedure  Decrement} 


Function  Missing(target,sen  :  integer; 

time.val  :  real)  :  boolean; 

begin 

Missing  :=  (val  <  missing_limit) 
or  (target  >  max.target) 
or  (sen  <>  sensor) 
or  (time  <  min.time) ; 

end;  {Function  Missing} 


Procedure  Echo _True_Data(target , sen  :  integer; 

time.val  :  real); 

begin 

if  (sen  =  sensor) 
and  (time  >=  min.time) 
and  (target  <=  max.target)  then 
begin 

while  (target  -  last.target)  >  1  do 
begin 

Write(trudata, ’  ’); 

Increment (last .target) ; 
end ;  {while} 

Increment (last .target) ; 
if  val  <  missing.limit  then 
begin 

Write(trudata,0:4) ; 
end  {if} 
else 
begin 

Write (trudata, target: 4) ; 
end;  {else} 
end;  {if} 

end;  {Procedure  Echo.True.Data} 

{***  Initializations  ♦*****************•*************************************}• 


Procedure  Init.Program;  {System  Specific} 

var 

i, j .k.atr.var.set  :  integer; 

filename  :  packed  array[1..12]  of  char; 


{Select  data  file  to  cluster} 

Readln(f ilename) ; 

{Determine  attribute  and  clock  biases} 
for  i  :=  1  to  attributes  do 
Readln(bias [i] ) ; 

Readln(tbias) ; 

{Determine  observation  sensor} 

sensor  :=  Ord(f ilename [2])  -  Ord(’O’); 

{Set  variances} 

var_set  :=  Ord(f ilename [3] )  -  Qrd(’O’); 
for  i  :=  1  to  attributes  do 
for  j  :=  1  to  attributes  do 
R_k[i,j]  :=  0.0; 

R_k[l ,1]  :=  10.0; 

R_k[2,2]  :=  1.0; 

R_k[3,3]  :=  Exp(-(var_set+2)*Ln(10.0)); 

R_k[4,4]  :=  Exp(-(var_set+2)*Ln(10.0)); 
for  atr  :=  1  to  attributes  do 
R[atr]  :=  Sqrt(R_k[atr,atr]); 
for  i  :=  1  to  3  do 

for  j  :=  1  to  nterms  do 
for  k  :=  1  to  nterms  do 
Rvar[i, j ,k]  :=  0.0; 
for  atr  :=  1  to  attributes  do 
begin 

Rvar [l , atr , atr]  : =  R_k [atr , atr] ; 

Rvar  [2 ,  atr ,  atr]  :  =  R.k  [atr ,  atr]  ; 

Rvar[2,atr+attributes,atr+attributes]  :=  R_k [atr , atr] ; 
end;  {for  atr} 

Rvar [3, 1 , 1]  :=  R_k[l,l]; 

Rvar [3, 2, 2]  :=  R_k[2,2]; 
for  i  :=  3  to  5  do 
begin 

Rvar[3,i,i]  :=  R_k[3,3]; 

Rvar [3 , i+3 , i+3]  :=  R_k[4,4]; 
end;  {for  i} 

case  filename[S]  of  {Determine  number  r'*  targets} 


*v>, 

’V* 

:  max.target 

:=  5 

’x\ 

’X’ 

:  max.target 

:=  10 

’T’ , 

’t> 

:  max.target 

:=  20 

’L\ 

*1’ 

:  max.target 

:=  50 

’C*. 

*  c 1 

:  max.target 

:=  100 

end; 

{case} 

case  filename[6]  of  {Determine  time  span} 
’A’, ’a’  :  begin 

min.time  :=  50.0; 
max.time  :=  200.0; 
end;  {Case  A} 


■tfSe 

l"  L'- 


i  *1* . .. . 


’B’,’b’  :  begin 

min.time  :=  0.0; 

max.time  :=  200.0; 

end;  {Case  B> 

'C’.’c'  :  begin 

min.time  :=  0.0; 

max.time  :=  300.0; 

end;  {Case  C> 

'D’.’d'  :  begin 

min_time  :=  0.0; 

max.time  :=  100.0; 

end;  {Case  D> 

'E'.’e'  :  begin 

min.time  :=  0.0; 

max.time  :=  2500.0; 
end;  {Case  E> 

'F’.’f'  :  begin 

min_time  :  =  0.0; 

max.time  :=  360.0; 

end;  {Case  F> 
end;  {case} 

case  filename [7]  of  {Determine  missing  data  rate} 

'■’,’n'  :  missing_limit  :=  0.00; 

'Y’.’y'  :  missing_limit  :=  0.05; 

’X’.’x’  :  mis8ing_limit  :=  0.10; 

’Z’.’z’  :  missing_limit  :=  0.20; 

end;  {case} 

{Initialize  input  data  file} 

Connect ( tmdata , ’TMDATA  ’); 

Reset (tmdata) ; 

{Initialize  sensor  input  data  file} 

Connect(ssdata, 'SSDATA  ’); 

Reset (ssdata) ; 

{Initialize  target  input  data  file} 

Connect (tsdata, 'TSDATA  ’); 

Reset (tsdata) ; 

{Initialize  cluster  progress  output  file} 

Connect (tcldata, ’TCLDATA  ’ ) ; 

Rewrite(tcldata) ; 

Writeln(tcldata, 'Clusters  for  attributes: 

’Reinge,  Reinge  Rate,  Azimuth,  and  Elevation’); 

Writeln(tcldata) ; 

{Initialize  true  observation  output  file} 

Connect (trudata, 'TRUDATA  ’); 

Revrite(trudata) ; 

Writeln(trudata, 'Clusters  for  attributes:  ’, 

'Range,  Range  Rate,  Azimuth,  and  Elevation  (True)’); 
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{Initialize  cluster  state  output  file} 

Connect (csdata, 'CSDATA  ’); 

Revrite(csdata) ; 

end;  {Procedure  Init .Program} 

{***  State  Derivative  Functions/Procedures  **********************************} 


real; 

stacked.vector; 
stacked.vector) ; 


integer; 

real; 

matrix; 


xr  :=  FX[l]*ri; 

yr  : 

=  FX  [2]  *ri ; 

zr 

x3x  :=  -3*Sqr(xr) ; 

y3y 

:=  -3*Sqr(yr); 

z3z 

x3y  :=  -3*xr*yr; 

x3z 

:=  -3*xr*zr; 

y3z 

fvr[l,l]  :=  mf *(i 

!•  x3x) 

» 

Procedure  Deriv(tm 
FX 

var  DX 

type 

matrix  =  array  [1. .block, 1. .block]  of  real; 
var 

i,j, split 

r , ri , r2i ,r3i ,mf , xr , yr , zr , x3x , x3y , x3z , y3y , y3z , z3z 
fvr 
begin 

{General  factors} 

r  :=  Sqrt(Sqr(FX[l])+Sqr(FX[2] )+Sqr(FX[3] )) ; 
ri  :=  1/r;  r2i  :=  Sqr(ri);  mf  :=  -mu*r2i/r; 

{State  transition  factors} 

:=  FX  [3]  *ri; 

:=  -3*Sqr(zr); 
:=  -3*yr*zr ; 

:=  mf*(i  !  x3x) ; 
fvr[l,2]  :=  mf*x3y; 
fvr Cl, 3]  :=  mf*x3z; 
fvr [2,1]  :=  mf*x3y; 
fvr [2,2]  :=  mf*(l  +  y3y); 
fvr [2, 3]  :=  mf*y3z; 
fvr [3,1]  :=  mf*x3z; 
fvr[3,2]  :=  mf*y3z; 
fvr[3,3]  :=  mf*(l  +  z3z); 

{State  derivatives} 
for  i  :=  1  to  3  do 
begin 

DX[i]  :=  FX[i+3] ; 

DX[i+3]  :=  mf *FX[i] 
end;  {for  i} 

{State  transition  derivatives} 
split  :=  block*nest; 
for  i  :=  1  to  split  do 

DX[nest+i]  :=  FX[nest+split+i] ; 
for  i  :=  1  to  block  do 
for  j  :=  1  to  nest  do 
DX[split+i*nest+j] 


{Position  derivatives} 
{Velocity  derivatives} 


end;  {Procedure  Deriv} 


=  fvr[i,l]*FX[  nest+j] 

+  fvr[i,2]*FX[2*nest+j] 

+  fvr[i,3]*FX[3*nest+j]  ; 


i 
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{***  Integration  Functions/Procedures  ***************************************} 
Procedure  Initialize_RK78; 


rkf [1] 

= 

2/27; 

rkf  [2] 

1/9; 

rkf  [3] 

= 

1/36; 

rkf  [4] 

= 

1/12; 

rkf  [5] 

= 

1/6; 

rkf  [6] 

= 

1/24; 

rkf [7] 

= 

1/8; 

rkf [8] 

= 

5/12; 

rkf  [9] 

= 

-25/16; 

rkf  [10] 

= 

1/2; 

rkf [11] 

= 

1/20; 

rkf  [12] 

= 

1/5; 

rkf  [13] 

= 

1/4; 

rkf [14] 

= 

5/6; 

rkf [15] 

= 

-25/108; 

rkf  [16] 

= 

125/108; 

rkf [17] 

= 

125/54; 

rkf [18] 

= 

-65/27 ; 

rkf  [19] 

= 

1/6; 

rkf [20] 

= 

13/900; 

rkf [21] 

= 

31/300; 

rkf  [22] 

= 

-2/9; 

rkf  [23] 

= 

61/225; 

rkf [24] 

= 

2/3; 

rkf  [25] 

= 

67/90; 

rkf [26] 

= 

-53/6; 

rkf [27] 

= 

-107/9; 

rkf  [28] 

= 

704/46; 

rkf  [29] 

= 

1/3; 

rkf  [30] 

= 

-1/12; 

rkf  [31] 

= 

23/108; 

rkf [32] 

= 

-19/60; 

rkf [33] 

= 

-91/108; 

rkf  [34] 

= 

17/6; 

rkf  [35] 

= 

311/54; 

rkf  [36] 

= 

-976/135; 

rkf  [37] 

= 

45/164; 

rkf [38] 

= 

18/41; 

rkf  [39] 

= 

2133/4100; 

rkf  [40] 

= 

45/82; 

rkf  [41] 

= 

2383/4100; 

rkf  [42] 

= 

-341/164; 

rkf  [43] 

= 

-301/82; 

rkf  [44] 

= 

4496/1025; 

rkf [45] 

= 

3/205; 

rkf  [46] 

= 

-3/41 ; 

rkf  [47] 

= 

-6/41; 

rkf  [48] 

= 

33/164; 

rkf  [49] 

= 

12/41; 

rkf [50] 

= 

-1777/4100; 

rkf  [51] 

= 

2193/4100; 

rkf  [52] 

= 

51/82; 

rkf  [63] 

= 

-341/164; 

rkf [54] 

= 

-289/82; 

rkf  [56] 

= 

4496/1025; 

rkf [56] 

= 

9/280; 

rkf [57] 

= 

41/840; 

rkf  [58] 

= 

9/35; 

rkf [59] 

= 

34/105; 

end;  {Procedure  Initialize  JUC78} 

Procedure  RK78(var  x  :  8tacked_vector ; 
var  t,dt  :  real; 
tout  :  real; 
relerr.abserr  :  real; 

neqn  :  integer; 
var  if lag  :  integer) ; 
label  1,2,3; 
const 

bup  =  2 . 821 109907456E+12 ; 
bio  =  1.6815125390625E-11; 
var 

dtfix.dtfail  :  boolean; 

i.nrej ,nrejt,nstp  :  integer; 

dtold, delt, tO, rer, scale, ae.rte.te.xmag, pet  :  real; 

fO, fl, 12, 13, f4, f5, f6,f7,f8,f9,fl0,f 11,112, xO  :  stacked_vector ; 

begin 

nrejt  :=  0; 
nstp  :=  0; 

if  (abserr  =  0)  and  (relerr  =  0)  then  {Set  flag  if  fixed  step  node  desired} 
dtfix  :=  true 

else 

dtfix  :=  false; 
dtold  :=  dt; 


l:dtfail  :=  falsa; 
nre j  : =  0 ; 

{Reset  stap  siza  if  this  sill  put  t  greater  than  tout} 
delt  :=  tout  -  t; 
if  Abs(dt)  >=  Abs(dalt)  then 
dt  : =  delt 
else 

if  Abs(2*dt)  >=  Abs(delt)  then 
dt  :=  delt/2; 

if  Abs(dt)  >=  2.84E-14*Abs(t)  then 
begin 

{First  Evaluation} 
tO  :=  t; 

for  i  :=  1  to  neqn  do 
xO[i]  :=  x[i]  ; 

Deriv(t,x,fO) ; 

{Second  Evaluation} 

2:  t  :=  tO  +  rkf[l]*dt; 
for  i  :=  1  to  neqn  do 

x[i]  :=  rkf [l}*fO[i]*dt  +  xO[i]; 

Deriv(t,x,f 1) ; 

{Third  Evaluation} 

t  :=  tO  +  rkf [2] *dt ; 
for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf  [3]*f0[i]  +  rkf  [4]*fi[i])*dt  +  xO[i]; 

Deriv(t ,x,f 2) ; 

{Fourth  Evaluation} 

t  :=  tO  +  rkf  [B]  *dt ; 
for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [6]*f0[i]  +  rkf  [7]*f2[i] )*dt  +  xO[i]; 

Deriv(t,x,f3) ; 

{Fifth  Evaluation} 

t  :=  tO  +  rkf [8] *dt ; 
for  i  :=  1  to  neqn  do 

xCi]  :=  (rkf [8]*f0[i]  +rkf [9]*(f2[i]  -  f3[i]))*dt  +  xO[i]; 
Deriv(t,x,f4) ; 

{Sixth  Evaluation} 

t  :=  tO  +  rkf[10]+dt; 
for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [ll]*f0[i]  +  rkf [12]*f4[i]  +  rkf [13]*f3[i] )*dt  +  xO[i]; 
Deriv(t,x,f5) ; 

{Seventh  Evaluation} 

t  :=  tO  +  rkf[14]*dt; 
for  i  :=  1  to  neqn  do 

xCi]  :=  (rkf [I5]*f0[i]  +  rkf L18]*f3[i]  +  rkf [17]vfS[i] 

+  rkf[l8]*f4Ci])*dt  +  x0[i}; 

Deriv(t,x,f8) ; 

{Eighth  Evaluation} 

t  ;=  tO  +  rkf[19]*dt; 


I  aU.^  i  J  »«.- U‘rf.rM"*.*,,*.m*fil‘t.l^1l '  ..  ,1*5.4'  «_•’* 


rkf  [31]*f3[i] 
rkf  [35]*f5[i] 


rki [32] *f6[i]  +  rkf  [33]*f0[i] 
rkf  [36]*f4[i])*dt  +  x0[i]; 


rkf[38]*f9[i]  +  rkf  [39]  *f  6  [i]  +  rkf  [40]*f7[i] 
rkf  [42]  *f  3  [i]  +  rkf  [43]*f5[i] 


for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [20]*f6[i]  +  rkf  [21]*f0[i]  +  rkf [22] *f5[i] 

+  rkf [23]*f4[i])*dt  +  xO[i]; 

Deriv(t,x,f7) ; 

{linth  Evaluation} 

t  :=  tO  +  rkf[24]*dt; 
for  i  :=  1  to  neqn  do 

xCi]  :=  (rkf  [26]*f6[i]  +  2*f0[i]  +  3*f7[i]  +  rkf  [26]  *f 3  [i] 

+  rkf [27]*f5[i]  +  rkf [28] *1 4  [i] ) *dt  +  xO[i]; 
Deriv(t,x,f8) ; 

{Tenth  Evaluation} 

t  :=  tO  +  rkf [29] *dt ; 
for  i  :=  1  to  neqn  do 

xti]  :=  (rkf [30]*f8[i]  +  rkf  [3i]*f3[i]  +  rkf  [32]*f6[i]  +  rkf[ 
+  rkf  [34]  *f  7  [i]  +  rkf  [35]*f5[i]  +  rkf [36] *f4[i] )*dt  + 
Deriv(t,x,f9) ; 

{Eleventh  Evaluation} 
t  :=  tO  +  dt; 
for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [37]*f8[i]  +  rkf [38]*f9[i]  +  rkf [39] *f 6 [i]  +  rkf[ 
+  rkf [41]*f0  [i]  +  rkf  [42]  *f  3  [i]  +  rkf  [43]  *f  5  [i] 

+  rkf  [44]  *f4  [i]  )  *dt  +  xO[i]  ; 

Deriv(t,x,f 10) ; 

{Twelfth  Evaluation} 
t  :=  tO; 

for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [45]*(f0[i]  -  f6[i])  +  rkf [46]* (f 7 [i]  -  f8[i]) 

+  rkf [47]* (f 5 [i]  -  f9[i]))*dt  +  xO[i]; 

Deriv(t,x,fll) ; 

{Thirteenth  Evaluation} 
t  :=  tO  +  dt; 
for  i  :=  1  to  neqn  do 

x[i]  :  =  (rkf [48]*f8[i]  +  rkf [49]*f9[i]  +  rkf [60]*f0[i]  +  rkf[ 
+  rkf [62]*f7[i]  +  f 11 [i]  +  rkf [53]*f3[i]  +  rkf[ 

+  rkf [56]*f4[i] )*dt  +  x0[i]; 

Deriv(t,x,f 12) ; 

{Compute  state  at  t+dt} 
for  i  :=  1  to  neqn  do 

x[i]  :=  (rkf [56]* (f 8 [i]  +  f9[i])  +  rkf [57]* (f 11 [i]  +  fl2[i]) 

+  rkf [58]*(f6[i]  +  f7[i])  +  rkf [59]*f5[i])*dt  +  x0[i]; 
if  not  dtfix  then 

begin  {Compute  max  local  truncation  error} 
rer  :=  RMax(relerr,2.572E-13) ; 
scale  :=  2/rer; 

ae  :=  scale  *  abserr  +  1.0E-14; 
rte  :=  0; 

for  i  :=  1  to  neqn  do 
begin 

te  :=  Abs(rkf [67] *(f0[i]  +  fl0[i]  -  fll[i]  -  f 12 [i] ) ) ; 


rkf [51]*f6[i] 
rkf  [64]  *f5[i] 


I 

I 


xnag  :=  Abs(x[i])  +  Abs(xO[i3) 
rte  :=  RMax(rte,te/xmag) ; 
end;  {for  i> 
rte  : =  rte  *  scale; 
if  rte  >=  1  then 

begin  {Reject  this  step> 
dtfail  :=  true; 
nrej  :=  nrej  +  1; 
nrejt  :=  nrejt  +  1; 
if  nrej  >=  10  then 
begin 

for  i  :=  1  to  neqn  do 
x[i]  xO  [i]  ; 

t  :=  tO; 

7; 


+  ae; 


if lag  := 
Goto  3; 
end  {if} 
else 
begin 
pet  :=  0 
if  rte  < 
pet  :  = 
dt  :=  pet 


025; 

bup  then 

0.9/Sqrt(Sqrt(Sqrt(rte))) ; 
dt; 


dt; 


dtold  :  = 

Goto  2; 
end;  {else} 
end;  {if} 

end;  {if  not  dtfix} 

{This  step  is  acceptable  -  eighth  order  evaluation} 
t  :=  tO  +  dt; 
nstp  :=  nstp  +  1; 
if  Abs(tout-t)  <=  1.0E-14  then 
begin 

dt  : =  dtold; 
if lag  :=  2; 

Goto  3; 
end; 

if  dtfix  then 
Goto  1; 
pet  :=  20; 
if  rte  >  bio  then 

pet  :=  0.9/Sqrt(Sqrt(Sqrt(rte))> ; 
if  dtfail  then 

pet  :=  RMin(pct.l.O); 
dt  :=  dt  *  pet; 
dtold  :=  dt; 

Goto  1; 
end;  {if} 
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{Check  for  too  small  a  step  size} 
if  Abs(delt)  <=  Abs(dt)  then 
begin  {Extrapolate} 

Deriv(t,x,fO) ; 
for  i  :=  1  to  neqn  do 
x[i]  :=  fO[i]*dt  +  x[i] ; 
t  :=  t  +  dt; 
dt  :=  dtold; 
if lag  :=  2; 
end 

else  {Attempted  to  use  too  small  a  step  size} 
if lag  :=  8; 

3: end;  {Procedure  RK78} 

{***  Integrator  Interface  Procedures  t***************************************} 

Procedure  Stack_Vector(x_vector  :  state.vector ; 

Phi_matrix  :  P_matrix; 
var  x_stack  :  stacked_vector) ; 

var 

i,j  :  integer; 
begin 

for  i  :=  1  to  nest  do 

x_stack[i]  :=  x_vector[i] ; 
for  i  :=  1  to  nest  do 
for  j  :*  1  to  nest  do 

x_stack[nest*i+j]  :=  Phi_matrix[i, j] ; 
end;  {Procedure  Stack .Vector} 

Procedure  Unstack_Vector(var  x.vector  :  state.vector; 

var  Phi.matrix  :  P_matrix; 

x.stack  :  stacked.vector) ; 

var 

i. j  :  integer; 
begin 

for  i  :=  1  to  nest  do 

x_vector[i]  :=  x_stack[i]; 
for  i  :=  1  to  nest  do 
for  j  :=  1  to  nest  do 

Phi_matrix[i, j]  :=  x_stack[nest*i+j] ; 
end;  {Procedure  Unstack.Vector} 

Procedure  Get .Sensor; 
var 

time, result  :  integer; 

step  :  real; 

t  :  array  [1..2]  of  real; 

SSV  :  stacked.vector ; 

Phi  :  P.matrix; 
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begin 

for  time  :=  framel  to  -1  do 
Sxs[time]  :=  Sxs[time+l]; 
repeat 

Read(sedata, current .sensor) ; 
until  (current .sensor .time  =  f rame.time [0] ) 
and  (current.sensor .number  =  sensor); 

Sxs[0]  :=  current.sensor .values; 

Phi  :=  Identity; 

Stack_Vector(Sxs [0] ,Phi,SSV); 
t [1]  :=  0.0;  t[2]  :=  tbias;  step  :=  tbias; 

RK78(SSV,t [1] ,8tep,t[2] .small, zero, nstack, result) ; 
Unstack_Vector(Sxs [0] ,Phi,SSV) ; 
end;  {Procedure  Get.Sensor} 


Procedure  Get .Observations ; 


t ime , nt ime , obs , atr 


integer; 

real; 


begin 

{Shift  previous  data} 

for  time  :=  framel  to  0  do 
begin 

ntime  :=  time  +  1; 

atr .limits [time]  :=  atr.limits Cntime] ; 
as signed [time]  : =  as signed [ntime] ; 
end;  {for  time} 
for  time  :=  framel  to  -1  do 
begin 

ntime  :=  time  +  1; 
nr.obsCtime]  :=  nr.obs [ntime] ; 
nr_unas8igned[time]  :=  nr.unassigned [ntime] ; 
nr .clusters [time]  :=  nr .clusters [ntime] ; 
nr_active[time]  :=  nr.active [ntime] ; 
nr_inactive[time]  :=  nr.inactive [ntime] ; 
attr[time]  :=  attr[ntime]; 
status [time]  :=  status [ntime] ; 
targets[time]  :=  targets [ntime] ; 
observation[time]  :=  observation [ntime] ; 
frame.t ime [time]  :=  frame.time [ntime] ; 
end;  {for  time} 

{Read  new  observations} 
obs  : =  1 ; 

attr[0,obs]  :=  next; 
targets [0, obs]  :=  next.target; 
frame.time[0]  :=  attr[0,obs,0] ; 
repeat 

E0I  :=  EOF(tmdata); 

Echo _True_Data(target8[0, obs] ,sen.nr,attr[0,obs,0j .missing.flag) ; 


w 


r 


ww 


to 

to 

as 


il  not  Missing(targets[0,obs] ,sen_nr ,attr[0,ob8,0] .missing.llag)  then 
begin 

lor  atr  :=  1  to  attributes  do 
begin 

atr.limits [0, 1 ,atr]  :=  RMin(atr.limits[0, 1 ,atr] , attr [0, obs, atr] ) ; 
atr .limits [0,2, atr]  :=  RMax(atr_limits[0,2,atr] ,attr [0,obs,atr] ) ; 
end;  {lor  atr} 

Increment (obs) ; 
end;  {il  not  Missing} 
il  not  E0I  then 
begin 

Read(tmdata, current _obs) ; 
target 8 [0, obs]  :=  current _obs .target; 
sen.nr  :=  current _obs . sensor; 
attr[0,obs,0]  :=  current _obs .time; 
lor  atr  :=  1  to  attributes  do 

attr [0 , obs , atr]  :=  current _obs . obs [atr] 

+  (bias [atr]  +  current _obs . error [atr] )*R [atr] ; 
missing_llag  :=  current _obs .missing; 
end;  {il  not  E0I} 

until  ( attr [0, obs, 0]  >  lrame_time[0])  or  E0I; 
last .target  :=  0; 

E0I  :=  E0I  or  ( attr [0, obs, 0]  >  max.time); 

Uriteln(trudata) ; 
il  not  E0I  then 
begin 

Write(trudata.attr[0,obs,0] :7: 1) ; 
next  :=  attr[0,obs]; 
next .target  : =  targets [0 , obs] ; 
end;  {il  not  E0I} 
nr_obs[0]  : =  obs  -  1; 
nr.unassigned[0]  :=  nr_obs[0]; 
nr_inactive[0]  :=  nr_inactive[-l] ; 
nr .active [0]  :=  nr .active [-1] ; 
nr .clusters [0]  :=  nr.active[0]  +  nr_inactive[0] ; 
end;  {Procedure  Get .Observations} 


Procedure  Initialize.Clustering; 


missed  :  boolean; 

time, ntime, cluster, obs, atr  :  integer; 
noise.ltime  :  real; 

begin 

{Initialize  attribute  minimums,  maximuas} 
lor  atr  :=  1  to  attributes  do 
begin 

atr_limits[l , 1 ,atr]  :=  big; 
atr .limits [1,2, atr]  :=  -big; 
end;  {lor  atr} 
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lor  obs  :=  1  to  clusters  do 
assignedtl ,obs]  :=  false; 
for  tine  : =  0  downto  franel  do 
begin 

ntine  :  =  tine  +  1; 

atr_linits[tine]  :=  atr.limits [ntine] ; 
as signed [tine]  :=  as signed [ntine] ; 
end;  {for  time} 

{Initialize  cluster  status} 
for  tine  :  =  frame 1  to  0  do 
begin 

nr .clusters [tine]  :=  0; 
nr .unassigned [time]  :=  0; 
nr.obs [time]  :=  0; 
nr_active[time]  :=  0; 
nr .inact  ive  [t  ime]  :  =  0 ; 

targets  [time, 0]  :=  0;  {*  For  missing  observations  *} 

frame_time[time]  :=  -1.0;  {*  For  output  buffering  *} 

for  cluster  :=  1  to  clusters  do 
begin 

observation[time, cluster]  :=  0; 
status [tine, cluster]  :=  0; 
end;  {for  cluster} 
end;  {for  time} 

{Find  first  measurement} 

It ime  :=  -1.0; 
last.target  :=  0; 
repeat 

ReadCtmdata, current _obs) ; 
next .target  :=  current _obs. target; 
sen.nr  :=  current _obs . sensor; 
next [0]  : =  current.obs . time ; 
for  atr  :=  1  to  attributes  do 

next[atr]  :=  current.obs. obs  [atr] 

+  (bias[atr]  +  current.obs . error [atr] )*R [atr] ; 
missing.flag  :=  current.obs. mis sing; 

if  next [0]  >  ltime  then  {*  Output  true  data  *} 

begin 

ltime  :=  next[0]; 

Uriteln(trudata) ; 

Write(trudata,ltime:7: 1) ; 
end;  {if} 

nissed  :=  Missing(next_target , sen.nr .next [0] .missing.flag) ; 
if  nissed  then 

Echo_True_Data(next .target , sen.nr , ltime .missing.flag) ; 
until  not  missed; 


for  time  :=  1  to  2  do 

begin  {Got  minimum  data  frames  minus  1  to  start  cluster} 
Get .Observations ;  {Read  remainder  of  observation  frame} 
Get .Sensor;  {Read  sensor  state} 

end;  {for  time} 

repeat  {Read  beginning  of  true  target  state  data  frame} 
Read(tsdata,current_target) ; 
until  (current .target. time  =  frame_time[0] ) ; 
end;  {Procedure  Initialize.Clustering} 


{***  Major  Procedures  *******************************************************} 


Procedure  Input .Data; 


time,obs,atr  :  integer; 
begin 

Start _Timer(l) ;  {Timing} 

Get .Observations ; 

Get .Sensor ; 

{Calculate  span  minimums  and  maximums  and  scale  factors} 
span.limits  :=  atr.limits [frame 1] ; 
for  atr  :=  1  to  attributes  do 
for  time  :=  framel+1  to  0  do 
begin 

span.limits [1 , atr]  :=  RMin(span_limits [l.atr] .atr.limits [time , 1 .atr] ) ; 
span.limits [2, atr]  :=  RMax(span_limits [2, atr] .atr.limits [time, 2, atr] ) ; 
end;  {for  time} 

Stop.Timer(l) ;  {Timing} 

end;  {Procedure  Input .Data} 


{***  Estimation  Initialization  Procedures  ***********************************} 


Procedure  Initialize.Estimation; 


i.j  :  integer; 
begin 

{Identity  matrix} 

for  i  :=  1  to  nest  do 
begin 

for  j  :=  1  to  nest  do 
Identity [i, j]  :=  0.0; 

Identity [i.i]  :=  1.0; 
end;  {for  i} 

{State  covariance} 

Q.k  :=  Identity; 
for  i  :=  1  to  nest  do 
t|_k[i,i]  :=  1.0; 
step  :=  10.0; 

end;  {Procedure  Initialize.Estimation} 
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{***  Estimation  Propagation  Procedures  **************************************} 

Function  ArcTan2(num,den  :  real)  :  real; 
var 

ansver  :  real; 
begin 

answer  :=  ArcTan(num/den) ; 
if  den  <  0  then 

answer  :=  answer  +  pi; 
if  answer  >  pi  then 

answer  :=  answer  -  2.0*pi; 

ArcTan2  :=  answer; 
end;  {Function  ArcTan2} 

Procedure  Map_State(Xr  :  state.vector ; 

var  Ov  :  obs_vector); 

begin 

Ov[l]  :=  Sqrt(Sqr(Xr [1] )+Sqr(Xr [2] )+Sqr(Xr [3] )) ; 

Ov  [2]  :=  (Xr  [l]*Xr[4]  +  Xrt2]*Xr[5]  +  Xr  [3]  *Xr  [6]  )/0v[l]  ; 

0v[3]  :=  ArcTan2(Xr[2] ,Xr[l]); 

Ov [4]  :=  ArcTan(Xr[3]/Sqrt(Sqr(Xr[l])+Sqr(Xr[2] ))) ; 
end;  {Procedure  Map_State> 

Procedure  Calculate_H(var  H  :  H_matrix; 

Xr  :  s>,ate_ vector; 
var  Ov  :  obs.vector) ; 

var 

i  :  integer; 

rhoi,rhoi2,varrho,varrhoi,varrhoi2,rhoi2_vari  :  real; 
begin 

{Calculate  coefficients} 

Map_State(Xr,Ov)  ; 
rhoi  :=  l/Ov[l]; 
rhoi2  :=  Sqr(rhoi); 

varrho  :=  Sqrt(Sqr(Xr [1] )  +  Sqr(Xr[2])); 
varrhoi  :=  1/ varrho; 
varrhoi2  :=  Sqr ( varrhoi) ; 
rhoi2_vari  :=  rhoi2*varrhoi; 

{Form  H  matrix} 

for  i  :=  i  to  block  do 
begin 

H[l,i]  :=  Xr[i]*rhoi; 

H[l,i+block]  :=  0.0; 

H [2 , i]  :=  (Xr[i+block] *0v[l]  -  Xr[i]*0v[2] )*rhoi2; 

H{2,i+block]  :=  H[l,i); 

H[3,i+block]  :=  0.0; 

H[4,i+block]  :=  0.0; 
end;  {for  i} 

H[3,l]  :=  -Xr [2] *varrhoi2; 
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H[3,2]  :=  Xr[l]*varrhoi2; 

H [3 , 3]  :=  0.0; 

H[4,l]  :=  -Xr[l]*Xr[3]*rhoi2_vari; 

H[4,2]  :=  -Xr [2] *Xr[3] *rhoi2_vari; 

H[4,3]  :=  varrho*rhoi2; 
end;  {Procedure  Calculate.!!} 

Procedure  Hap_State_to_Attribute_Space(Xv  :  state.vector ; 

Pm  :  P .matrix; 
var  0v,0P  :  obs.vector), 

var 

i,j,k  :  integer; 

Xr,SP  :  state.vector ; 

H,H_P  :  H .matrix; 

M  :  R .matrix; 
begin 

for  j  :=  1  to  nest  do 

Xr[j]  :=  Xv[j]  -  Sxs [0 , j] ; 

Calculate_H(H,Xr ,0v) ; 
for  i  :=  1  to  attributes  do 
for  j  :=  1  to  nest  do 
begin 

H_P[i, j]  :=  0.0; 
for  k  :=  1  to  nest  do 

H_P[i, j]  :=  H_P[i, j]  +  H[i,k]*Pm[k, j] ; 
end;  {for  j} 

M  :=  R_k; 

for  i  :=  1  to  attributes  do 
for  j  :=  1  to  attributes  do 
for  k  :=  1  to  nest  do 

M[i,j]  :=  M[i , j]  +  H_P[i,k]*H[j ,k] ; 
for  i  :=  1  to  attributes  do 
0P[i]  :=  3.0*Sqrt(M[i,i] ) ; 
end;  {Procedure  Map_State_to_Attribute> 

Function  ArcSin(arg  :  real)  :  real; 
begin 

ArcSin  :=  ArcTan(arg/Sqrt(l .O-Sqr(arg))) ; 
end;  {Function  ArcSin} 

Function  ArcCos(arg  :  real)  :  real; 
begin 

ArcCos  :=  pi/2.0  -  ArcTan(arg/Sqrt(l . 0-Sqr (arg) ) ) ; 
end;  {Function  ArcCos} 

Procedure  Forecast ; 
var 

cluster, i,j ,k, result  :  integer; 
t  :  array  [1..2]  of  real; 
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stacked.vector ; 
P_matrix; 


SSV 

Phi,Phi_P 
begin 

Start_Timer(2) ;  {Timing} 

lor  cluster  :=  1  to  nr .clusters [-1]  do 
il  status [0 , cluster]  <  0  then 

begin  {Integrate  state  and  state  transition  matrix} 

Phi  :=  Identity; 

Stack_Vector(Sx [cluster] , Phi, SSV) ; 

t [13  :=  0.0; 

t[2]  :=  frame_time[0]  -  lrame_time [status [0, cluster]] ; 
status [0, cluster]  :=  -1; 

RK78(SSV,t [1] ,step,t[2] , small, zero, nstack, result) ; 
Unstack_Vector(Sx[cluster] , Phi, SSV) ; 

{Propagate  state  covariance  matrix} 
lor  i  :=  1  to  nest  do 
lor  j  :=  1  to  nest  do 
begin 

Phi_P[i, j]  :=  0.0; 
lor  k  :=  1  to  nest  do 

Phi_P[i,j]  :=  Phi_P[i,j]  +  Phi[i,k]*P[cluster,k, j] ; 
end;  {lor  j} 

P [cluster]  :=  q_k; 
lor  i  :=  1  to  nest  do 
lor  j  :=  1  to  nest  do 
lor  k  :=  i  to  nest  do 

P [cluster, i,j]  :=  P[cluster,i, j]  +  Phi_P[i,k]*Phi[j ,k] ; 
Hap_State_to_Attribute_Space(Sx[cluster] ,P[cluster] , 

est [cluster] .gate [cluster]) ; 

end;  {il} 

Stop_Timer(2) ;  {Timing} 

end;  {Procedure  Forecast} 

Procedure  Calculate.Hetrics; 
label  i; 
var 

cluster ,obs ,atr ,mro«,missed_gates  :  integer; 

delta  :  real; 

lactor  :  obs_vector; 

nr_leasible  :  array  [1..2]  ol  ivector; 

begin 

Start .Timer (3) ;  {Timing} 

mrow  :=  0; 

{Calculate  standardization  lactors} 
lor  atr  :=  1  to  attributes  do 
begin 

lactor[atr]  :=  span.limits [2 ,atr]  -  span.limits [1 ,atr] ; 
il  Abs (lactor [atr])  <  zero  then 
lactor [atr]  :=  1.0 
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factor Catr]  : =  1 . O/factor Catr] ; 
end;  {for  atr} 

{Initialize  pre-assignment  indices} 
for  cluster  :=  1  to  clusters  do 
nr .feasible  Cl, cluster]  :=  0; 
nr_feasible[2]  :=  nr_f easibleCl] ; 
col_basis  :=  nr _f easibleCl] ; 
row.basis  :=  nr _f easibleCl] ; 

{Calculate  metrics} 

for  cluster  :=  1  to  nr .clusters C-l]  do 
if  status  CO, cluster]  =  -1  then 
begin 

Increment (mrow) ; 
convert Cmrow]  :=  cluster; 
for  obs  :=  1  to  nr.obsCO]  do 
begin 

metric [mrow, obs]  :=  0.0; 
mis8ed_gates  :=  0; 
for  atr  :=  1  to  attributes  do 
begin 

delta  :=  Abs(est Ccluster ,atr]  -  attr CO , obs , atr] ) ; 
if  delta  >  gate Ccluster, atr]  then 
if  (delta  >  2. 0*gate Ccluster, atr] ) 
or  (mi88ed_gates  >=  max.missed)  then 
begin 

metric  Cmrow , obs]  : =  bad ; 
goto  1; 

end  {if  missed} 
else 

Increment (missed .gates) ; 

metric Cmrow , obs]  :=  metricCm^ow.obs]  +  Sqr (factor Catr] *delta) ; 
end;  {for  atr} 

1:  if  metric Cmrow, obs]  <  bad  then 

begin 

Increment (nr .feasible  Cl ,mrow] ) ; 

Increment (nr_feasible[2, obs] ) ; 
col.basis Cmrow]  :=  obs; 
row.basis Cobs]  :=  mrow; 
end;  {if} 
end;  {for  obs} 
end;  {if  status} 

{Check  pre-assignment} 

for  cluster  :=  1  to  mrow  do 

if  (nr.feasible Cl .cluster]  >  1)  or 

(nr _feasibleC2, col.basis Ccluster]]  >  1)  then 
col.basis Ccluster]  :=  0; 
for  obs  :=  1  to  nr.obsCO]  do 
if  (nr_feasibleC2,obs]  >  1)  or 
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(nr_feasible[2,row_basi8[obs]]  >  1)  then 
row.basis [obs]  : =  0 ; 

{Assign  dummy  values,  as  necessary} 

ior  cluster  :=  nr_clusters[-l]+l  to  nr_obs[0]+nr_inactive[-l]  do 
begin 

Increment (mrow) ; 
convert [mrow]  : =  cluster; 

status [0, cluster]  :=  0;  {Du 

for  obs  :=  1  to  nr_obs[0]  do 
metricOnrow.obs]  :=  bad; 
end;  {for  cluster} 

for  obs  :=  nr_obs[0]+l  to  nr_active[-l]  do 
begin 

targets [0, obs]  :=  -1;  {Dummy- 

for  cluster  :=  1  to  mrow  do 
metric [cluster, obs]  :=  bad; 
end;  {for  obs} 

Stop_Timer(3) ; 

end;  {Procedure  Calculate_Metrics} 


{Dummy  cluster} 


{Dummy-obs  ervat ion} 


{Timing} 


{***  State  Update  Procedures  ************************************ ************y 

Procedure  Invert (M  :  R_matrix; 

var  Mlnv  :  R_matrix) ; 

label  1; 
var 

i, j ,k,l, irow.icol.il  :  integer; 

determ, pivot, hold, sum.t.ab, big  :  real; 

index  :  array [1 . .nest , 1 .. 3]  of  integer; 

Procedure  SwapCvar  a,b  :  real); 
var 

hold  :  real ; 
begin 

hold  :  =  a; 
a  :=  b; 
b  :=  hold; 

end;  {Procedure  Swap} 

{Gauss- Jordan  inversion} 
begin 

for  i  :=  1  to  attributes  do 
index  [i ,  3]  :  =  0 ; 
determ  :=  1; 

for  i  :=  1  to  attributes  do 

begin  {Search  for  largest  element} 
big  :=  0; 

for  j  :=  1  to  attributes  do 
begin 

if  index  [j  ,3]  <>  1  then 
begin 
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lor  k  :=  1  to  attributes  do 


begin 

il  index  [k,  3]  >  1  then 
begin 

writeln(’ ERROR:  Matrix  singular'); 
goto  1; 
end; 

il  index [k, 3]  <  1  then 
il  Abs(M[j,k])  >  big  then 
begin 
iron  :=  j; 
icol  :=  k; 
big  :=  Abs(M[j ,k] ) ; 
end;  {il} 
end;  {lor  k> 
end;  {il} 
end;  {lor  j} 

Increment (index [icol , 3] ) ; 
index  [i ,  1]  ;  =  irow ; 
index [i, 2]  :=  icol; 

{Interchange  rows  to  put  pivot  on  diagonal} 
il  iron  <>  icol  then 
begin 

determ  :=  -determ; 
lor  1  :=  1  to  attributes  do 
Swap (M [irow, 1] ,M[icol,l]) ; 
end;  {il  iron  <>  icol} 

{Divide  pivot  row  by  pivot  column} 
pivot  :=  M [icol, icol] ; 
determ  :=  determ  *  pivot; 

M [icol , icol]  : *  1 ; 
lor  1  :=  1  to  attributes  do 

M[icol,l]  :=  M[icol,l]  /  pivot; 

{Reduce  nonpivot  rows} 

lor  11  :=  1  to  attributes  do 
begin 

il  11  <>  icol  then 
begin 

t  :=  M[ll , icol] ; 

M [11 .icol]  :=  0; 
lor  1  :=  1  to  attributes  do 

M[ll ,1]  :=  M[ll ,1]  -  M[icol,l]  *  t; 
end;  {il  11  <>  icol} 
end;  {1 or  11} 
end;  {lor  i} 

{Interchange  columns} 

1 or  i  :=  1  to  attributes  do 
begin 

1  :=  attributes  -  i  +  1; 


Loswwsxjw 


$ 

:s 


cuw cvoyovcyKJWV' 


if  index [1,1]  <>  index [1,2]  then 
begin 

irov  : =  index [1,1]; 
icol  : =  index [1,2]; 
for  k  :=  1  to  attributes  do 
svap(M[k,irov] ,M[k,icol] ) ; 
end;  {if  index} 
end;  {for  i} 

for  k  :=  1  to  attributes  do 
if  index Dc, 3]  <>  1  then 
begin 

sriteln( ’ ERROR :  Matrix  singular'); 

goto  1; 

end; 

for  i  :=  1  to  attributes  do 
for  j  :=  1  to  attributes  do 
MInv[i,j]  :=  M[i, j] ; 
l:end;  {Procedure  Invert} 

Procedure  Cholesky (M  :  P_matrix; 

var  S  :  P_matrix); 

var 

i.j.k  :  integer; 
sum  :  real; 
begin 

for  i  :=  1  to  nest  do 


begin 
for  j 


:=  1  to  i-1  do 


S[j,i]  :=  0.0; 
sum  : =  0.0; 
for  k  :=  1  to  i-1  do 

sum  :=  sum  +  Sqr(S[i,k]); 

S[i,i]  :=  Sqrt(M[i, i]  -  sum); 
for  j  :=  i+1  to  nest  do 
begin 

sum  :=  0.0; 

for  k  :=  1  to  i-1  do 

sum  :=  sum  +  S[i,k]*S[j ,k] ; 

S [ j , i]  :=  (M[i, j]  -  sum)/S[i,i] ; 
end;  {for  i} 
end;  {for  i} 

end;  {Procedure  Cholesky} 


Procedure  Update.Estimates ; 
var 

cluster , i , j , k , obs 
Oref , od 
dS.dSx 
H 


integer; 
obs_vector ; 
state_vector ; 
H_matrix; 
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K_k,Fbar,F_M  :  K .matrix; 

N,MI  :  R.matrix; 

Vfbar , What , I _FMF , Lambda , Phat  :  P_matrix; 
begin 

Start_Timer(6) ; 

lor  cluster  :=  1  to  nr .clusters [0]  do 
if  (status [0, cluster]  =  -1)  and 
(nr .missing [cluster]  =  0)  then 
begin 

{Calculate  relative  state  vector} 
lor  i  :=  1  to  nest  do 

dS[i]  :=  Sx [cluster, i]  -  Sxs[0,i]; 

{Calculate  H} 

Calculate_H(H,dS,Orel) ; 

{Calculate  W  bar} 

Cholesky(P [cluster] ,Vbar) ; 

{Calculate  F  bar} 

lor  i  :=  1  to  nest  do 

lor  j  :=  1  to  attributes  do 
begin 

Fbar[i,j]  :=  0.0; 

lor  k  :=  i  to  nest  do  {Lower  triangular  multiplication} 
Fbar [i,  j]  :=  Fbar[i,j]  +  Wbar  [k ,  i]  *H  [  j  ,k]  ; 
end;  {lor  j} 

{Calculate  H  matrix} 

MI  :=  R.k; 

lor  i  :=  1  to  attributes  do 
lor  j  :=  1  to  attributes  do 
lor  k  :=  1  to  nest  do 

MI[i,j]  :=  MI[i,j]  +  Fbar[k,i]*Fbar[k,j] ; 

Invert (MI, M) ; 

{Calculate  Lambda  matrix} 

lor  i  :=  1  to  nest  do 

lor  j  :=  1  to  attributes  do 
begin 

F_M[i, j]  :=  0.0; 

lor  k  :=  1  to  attributes  do 

F_M[i, j]  :=  F_M[i, j]  +  Fbar[i,k]*M[k,j]; 
end;  {lor  j} 

I..FHF  :=  Identity; 
lor  i  :=  1  to  nest  do 
lor  j  :=  i  to  nest  do 

lor  k  :=  1  to  attributes  do 

I_FMF[i,j]  :=  I_FMF[i,j]  -  F_M[i,k]*Fbar[j ,k] ; 
Cholesky( I _FMF, Lambda) ; 

{Calculate  V  hat} 

1 or  i  :=  1  to  nest  do 
begin 

lor  j  :=  1  to  i  do 


{Timing} 
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begin 

What[i,j]  :=  0.0; 
lor  k  :=  1  to  i  do 

What[i,j]  :=  What[i,j]  +  Vbar [i,k]*Lambda[k,j] ; 
end;  {lor  j} 
lor  j  :=  i+1  to  neat  do 
lfhat[i,j]  :=  0.0; 
end;  {lor  i> 

{Calculate  P  hat} 

lor  i  :=  1  to  nest  do 
lor  j  :=  1  to  i  do 
begin 

Phat[i,j]  :=  0.0; 

1 or  k  :=  1  to  i  do 

Phat[i,j]  :=  Phat[i,j]  +  What[i,k]*What[j ,k] ; 

PhatCj.i]  :=  Phat[i, j] ; 
end;  {lor  j} 

{Calculate  K  matrix} 

lor  i  :=  1  to  nest  do 

lor  j  :=  1  to  attributes  do 
begin 

K_k[i,j]  :=  0.0; 
lor  k  :=  1  to  i  do 

K_k[i,j]  :=  K_k[i,j]  +  Wbar[i,k]*F_M[k, j] ; 
end;  {lor  j} 

obs  :=  obs ervat ion [0, cluster] ; 
lor  i  :=  1  to  attributes  do 

od[i]  :=  attr [0 , obs , i]  -  OrelCi]; 
lor  i  :=  1  to  nest  do 
begin 

dSxCi]  :=  0.0; 

lor  j  :=  1  to  attributes  do 

dSx [i]  :=  dSxCi]  +  K_k[i, j]*od[j] ; 

Sx[cluster,i]  :=  SxCcluster , i]  +  dSxCi];  {Rectify  reference  state} 
end;  {for  i} 

P [cluster]  : =  Phat ; 
end;  {if} 

Stop_Timer(6) ;  {Timing} 

end;  {Procedure  Update.Estimates} 


Procedure  Assignment (n  :  integer; 

cost  :  M_matrix; 
var  ans  :  ivector) ; 

label  1,2; 
type 

rvector  =  array  [1. .clusters]  of  real; 

bvector  =  array  [1 . .clusters]  of  boolean; 

imatrix  =  array  [1. .clusters, 1. .clusters]  of  integer; 


{Hungarian  Method} 
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i.j.k.m  :  integer; 

delta  :  real; 

alpha, beta, slack  :  rvector; 

labels, exposed, nhbor.q  ;  ivector; 

labeled  :  bvector; 

mate  :  array  [1. 

A  ;  imatrix ; 

Procedure  Augment (v  :  integer); 
begin 

if  labels [v]  =  0  then 
begin 

mateCl.v]  :=  exposed [v] ; 
mate [2, exposed [v]]  :=  v; 
end  {if} 
else 
begin 

exposed[labels[v]]  :=  mateCl.v]; 
mate[l,v]  :=  exposed [v] ; 
mat e C2, expo sedCv]]  :=  v; 
augment (labels Cv] ) ; 
end;  {else} 

end;  {Procedure  Augment} 


:  integer; 

:  real; 

:  rvector; 

,Q  :  ivector; 

:  bvector; 

:  array  Cl. -2]  of  ivector; 
:  imatrix ; 
integer); 


Function  Modify 
label  1; 


boolean; 


i,j  :  integer; 

thetal,theta2  :  real; 
begin 

thetal  :=  big; 
for  j  ;=  1  to  n  do 
if  slackCj]  >  0  then 
thetal  :=  RMin(thetal, slack ( j] ) ; 
theta2  :=  thetal/2; 
for  i  :=  1  to  n  do 
if  labeled Ci]  then 
alpha Ci]  :=  alpha Ci]  +  theta2 
else 

alphaCi]  :=  alphaCi]  -  theta2; 
for  j  :=  1  to  n  do 
if  slackCj]  =  0  then 
betaCj]  :=  betaCj]  -  theta2 
else 

betaCj]  :=  betaCj]  +  theta2; 
for  j  :=  1  to  n  do 
if  slack Cj)  >  0  then 
begin 

slackCj]  :=  slackCj]  -  thetal; 
if  slackCj]  =  0  then 


il  mate [2, j]  =  0  then 
begin 

exposed [nhbortj]]  :=  j; 
augment (nhbor [j] ) ; 

Modify  :=  true; 
goto  1; 
end  {if} 
else 
begin 

labels [mate [2, j]]  :=  nhbor [j]; 
labeled [mate [2, j]]  :  =  true; 

Q[mate[2,j]]  :=  1; 

A[nhbor[j] ,mate[2, j]]  :=  1; 
end;  {else} 
end;  {if} 

Modify  : =  false; 

1:  end;  {Function  Modify} 

Function  q.not .empty (var  index  :  integer)  :  boolean ; 
begin 

index  :=  0; 
repeat 

Increment (index) ; 

until  (q [index]  =  1)  or  (index  =  n) ; 
if  (index  =  n)  and  (q [index]  =  0)  then 
q_not .empty  :=  false 
else 

q.not .empty  :=  true; 
end;  {Function  q.not.empty} 
begin 

for  i  :=  1  to  n  do 
begin 

mate[l,i]  :=  0; 
alpha [i]  :=  0; 
labels [i]  :=  0; 
end;  {for  i} 
for  j  :=  1  to  n  do 
begin 

mate[2,j]  :=  0; 
beta[j]  :=  big; 
for  i  :=  1  to  n  do 

beta[j]  :=  RMin(beta[j] ,cost[i,j]) ; 
end;  {for  j} 

{Repeat  for  n  stages} 
for  m  :=  1  to  n  do 
begin 

for  i  :=  1  to  n  do 
for  j  ;=  1  to  n  do 
A[i,j]  ;=  0; 
for  i  :=  1  to  n  do 
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sew 


«xposed[i]  :=  0; 
for  j  :=  1  to  n  do 
alack [j]  :=  big; 
for  i  :=  1  to  n  do 
for  j  :=  1  to  n  do 

if  cost[i,j]  -  alphaCi]  -  betaCj]  <  zero  then 
if  mate[2,j}  =  0  then 
begin 

exposed [i]  : =  j ; 
end  {if} 
else 
begin 

A[i,mate[2, j]]  r=  l; 
end;  {else} 

{Construct  auxiliary  graph} 
for  i  :=  1  to  n  do 
begin 

qCi]  :=  0; 

labeled[i]  :=  false; 
end;  {for  i} 
for  i  :=  1  to  n  do 
if  mateCl.i]  =  0  then 
begin 

if  exposed [i]  <>  0  then 
begin 

augment (i) ; 
goto  2; 
end;  {if} 

q[i]  :=  1; 

labels  [i]  :  =  0 ; 
labeled [i]  :=  true; 
for  k  :=  1  to  n  do 
begin 

delta  :=  costCi.k]  -  alphaCi]  -  beta[k] ; 
if  delta  <  zero  then 
delta  :=  0.0; 

if  (0  <=  delta)  and  (delta  <  slack [k] )  then 
begin 

slackCk]  :  =  delta; 
nhborCk]  :=  i; 
end;  {if} 
end;  {for  k} 
end; {if} 

1:  while  q_not_empty(i)  do 
begin 

q[i]  :=  0; 

if  exposed [i]  <>  0  then 
begin 

augment (i) ; 
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goto  2; 
end; 

lor  k  :=  1  to  a  do 
begin 

delta  :=  cost[i,k3  -  alphaCi]  -  beta[k]; 
if  delta  <  zero  then 
delta  :=  0.0; 

if  (0  <=  delta)  and  (delta  <  slack [k] )  then 
begin 

slack [k]  :=  delta; 
nhborCk]  :=  i; 
end;  {if} 
end;  {for  k} 
for  j  :=  1  to  n  do 

if  (A[i,j]  =  1)  and  not  labeledCj]  then 
begin 

labels  [j]  :=  i; 
labeled[j]  :=  true; 

QCj]  :=  1; 

if  exposed [j]  <>  0  then 
begin 

augment(j) ; 
goto  2; 
end;  {if} 

for  k  :=  1  to  n  do 
begin 

delta  :=  costCj.k]  -  alphaCj]  -  betaCk]; 
if  delta  <  zero  then 
delta  :=  0.0; 

if  (0  <=  delta)  and  (delta  <  slackCk])  then 
begin 

slack [k]  :  =  delta; 
nhbor [k]  : =  j ; 
end;  {if} 
end;  {for  k} 
end;  {if} 
end;  {while} 
if  not  Modify  then 
goto  i; 

2:  end;  {for  m} 
ana  :=  mate [2]; 
end;  {Procedure  Assignment} 

Procedure  Perf ona_Cluster_Assignment ; 
var 

cluster , obs , k , atr .size , 
asize,col_index,row_index  :  integer; 
brow,row_convert  :  ivector; 

done  :  boolean; 


9 

& 

V* 


•‘.  .t-.v  .-i I'l  i'i  .MYir.'iU'i'rtri'iri'tMUi.vtVi'i'UijNiiWi'iyti'filV 


li 


begin 

Start _Tiaer(4) ; 

size  :=  IMax(nr_active[-l] ,nr_obs[0]  ) ; 

{Perform  pre-assignment  packing} 
row.index  :=  0; 
asize  :=  0; 

lor  cluster  :=  1  to  size  do 
if  col_bas is [cluster]  =  0  then 
begin 

Increment (as ize) ; 

Increment (row.index) ; 
row.convert [row.index]  :=  cluster; 
col.index  :=  0; 
for  obs  :=  1  to  size  do 
if  row.basis [obs]  =  0  then 
begin 

Increment(col_index) ; 

ametric [row.index, col_index]  :=  metr ic  [cluster, obs] ; 
end;  {if} 
end;  {if} 
if  asize  >  0  then 
begin 

As8ignment(asize,ametric,brow) ; 
col.index  :=  0; 

for  obs  :=  1  to  size  do  {Unpack  solution} 
if  row.basis [obs]  =  0  then 
begin 

Increment(col_index) ; 

row.basis [obs]  :=  row.convert [brow [col.index]] ; 
end;  {if} 
end;  {if} 
brow  ;=  row.basis; 

{Check  for  invalid  assignments  and  assign  targets  to  clusters} 
for  obs  :=  1  to  size  do 
begin 

cluster  :=  convert [brow [obs]] ; 
if  metric [brow [obs] ,obs]  <  bad  then 
begin  {Good  assignment} 
observation[0,clu8ter]  :=  obs; 
assigned[Olobs]  :=  true; 

Deer ement(nr_un assign ed[0]) ; 
nr .missing [cluster]  :=  0; 
end  {if} 
else 

if  status [0, cluster]  =  -1  then 
begin  {Kissing  assignment} 

{Cluster  has  been  propagated  maximum  number  of  times} 
if  nr .missing [cluster]  >=  max.prop  then 
begin  {Terminate  old  cluster} 


{Timing} 
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status [0, cluster]  :=  0; 

Increment (nr .inactive [0] ) ; 
Decrement(nr.active[0] ) ; 
end  {if} 

else 

begin  {Propagate  old  cluster} 

Increment (nr .missing [cluster] ) ; 
observation^, cluster]  :=  0; 
end;  {else} 
end;  {else} 
end;  {for  obs} 

Stop_Timer(4) ; 

end;  {Procedure  Perform.Cluster .Assignment} 

Procedure  Calculate_SP_Estimates(timel,obsl  :  integer; 

var  est.gate  :  0 .matrix) ; 

type 

vector  =  array  [1.. block]  of  real; 
var 

i,j,time2  :  integer; 

A , acc , acc.n , acc.p ,acc_t2 , Arrs , 

Arrs_dot_urho,Arrsv,bl ,b2, 
beta.l , bet a_2 , beta.min , beta  jnax , 

Cos.az , Cos.el , Cos.beta , Cos.beta.l , 

Cos _beta_2 , d , dl , d2 , d3 , dmax , dmin , 
mu.f actorl ,mu_f actor2 ,rl_dot_rs2 , 
rng.min , rng.max , Sin.az , Sin.el , Sqr_t21 , 
t2 1 , theta , vmag2 , vsmag2 , vs .dot  _urho , 
rmag2,rmag3  :  real; 
rsmag2  :  array  [1..2]  of  real; 

urho.rvec  :  vector; 

begin 

{Calculate  unit  range  vector} 

Cos.az  :=  Cos(attr [timel ,obsl ,3] ) ; 

Cos.el  :=  Cos (attr [timel ,obsl ,4] ) ; 

Sin.az  :=  Sin(attr [timel ,obsl ,3] ) ; 

Sin.el  :=  Sin(attr [timel ,obsl ,4] ) ; 
urho[l]  :=  Cos_el*Cos_az; 
urho[2]  :=  Cos_el*Sin_az; 
urho[3]  :=  Sin.el; 

{Compute  magnitudes;  dot  products} 
rmag2  :=  0.0;  rsmag2[l]  :=  0.0; 

vsmag2  :=  0.0;  vs.dot.urho  :=  0.0; 
for  i  :=  1  to  block  do 
begin 

j  :=  i  +  block; 

rvec[i]  :=  attr [timel ,obsl , 1] *urho  [i]  +  Sxs [timel , i] ; 

rmag2  :=  rmag2  +  Sqr(rvec[i] ) ; 

rsmag2[l]  :=  rsmag2[l]  +  Sqr (Sxs [timel , i] ) ; 
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vsmag2  :=  vsmag2  +  Sqr(Sxs [timel , j] ) ; 
v8_dot_urho  :=  vs_dot_urho  +  Sxs [timel, j]*urho[i] ; 
end;  {for  i> 

vmag2  :=  2.0*(Max_Energy  +  mu/Sqrt (rmag2) ) ; 
mu_factorl  :=  -mu/(Sqrt(rmag2)*rmag2) ; 
mu_factor2  :=  -mu/(Sqrt(rsmag2[l] )*rsmag2[l] ) ; 
acc_p  :=  0.0;  acc_t2  :=  0.0; 
for  i  :=  1  to  block  do 
begin 

acc  :=  mu_factorl*rvec[i]  -  mu_f actor2*Sxs [timel , i] ; 
acc_t2  :=  acc_t2  +  Sqr(acc) ; 
acc_p  :=  acc_p  +  acc*urho[i]; 
end;  {for  i> 

acc_n  :=  Sqrt(acc_t2  -  Sqr(acc_p)); 

Cos_beta_l  :=  ( at tr [timel ,obsl,2]  +  vs_dot_urho)/Sqrt(vmag2) ; 
beta_l  :=  ArcCos(Cos_beta_l) ; 

dl  :=  Sqrt(vmag2  -  Sqr(attr [timel ,obsl ,2]  +  vs_dot_urho) ) ; 
d2  :=  Sqrt(vsmag2  -  Sqr(vs_dot_urho)) ; 
for  time2  :=  timel+1  to  0  do 
begin 

{Calculate  time  interval} 

t21  :=  frame_time[time2]  -  frame_time [timel] ; 

Sqr_t21  :=  0.5*Sqr(t21) ; 

{Compute  magnitudes;  dot  products} 
rsmag2 [2]  : =  0 . 0 ; 

rl_dot_rs2  :=  0.0; 
for  i  :=  1  to  block  do 
begin 

rsmag2 [2]  :=  rsmag2[2]  +  Sqr(Sxs [time2, i] ) ; 
rl_dot_rs2  :=  rl_dot_rs2  +  rvec[i]*Sxs[time2,i] ; 
end;  {for  i} 

A  :=  1.0  +  mu_factorl+Sqr_t21 ; 

{Calculate  accelerations;  more  dot  products} 

Arrs  :=  0.0;  Arrs_dot_urho  :=  0.0; 
for  i  :=  1  to  block  do 
begin 

Arrsv  :=  A*rvec[i]  -  Sxsrtime2,i] ; 

Arrs  :=  Arrs  +  Sqr (Arrsv); 

Arrs_dot_urho  :=  Arrs_dot_urho  +  Arrsv*urho[i] ; 
end;  {for  i} 

Arrs  :=  Sqrt(Arrs); 

{Compute  angle  beta} 

Cos_beta_2  :=  Arrs_dot_urho/Arrs ; 

Cos_beta  :=  RMin(Cos_beta_l ,Cos_beta_2) ; 
beta_2  :=  ArcCos(Cos_beta_2) ; 
beta_max  :=  beta_l  +  beta_2; 
beta.min  :=  Abs(beta_l  -  beta_2); 

{Calculate  range  estimate  and  gate} 

bl  :=  Sqr(A)*rmag2  +  Sqr(t21)*vmag2  -  2. 0*A*rl_dot_rs2  +  rsmag2[2]; 


b2  :=  2.0*t21*Arrs*Sqrt(vmag2); 

eat [time2, 1]  :=  Sqrt(bl  +  b2*Cos_beta_l) ; 

gate[tiae2, 1]  :=  Abs(Sqrt(bl  +  b2*Cos(beta_max))  -  est  [time2, 1]  )  ; 
gate[time2, 1]  :=  RHax(Abs(Sqrt(bl  +  b2*Cos(beta_Bin))  -  est [tiBe2, 1] ) , 

gate [time2 , 1] )  +  3.0*R[1]; 

{Calculate  range  rate  estiBate  and  gate} 
d3  :=  Sqr_t21*acc_n; 
dm? i  :=  t21*(dl+d2)  +  d3; 
dmm  :=  t21*Abs(dl-d2)  -  d3; 
d  :=  0.5*(dmax  +  dmin) ; 

est [tiae2,2]  :=  (Sqr(est[time2,l])  -  attr[timel,obsl,l] 

♦Sqrt (Sqr (est [time2 , 1] )  -  Sqr(d)))/(est[time2,l]*t21) ; 
rng_min  :=  est[time2,l]  -  gate [time2 , 1] ; 
mg_max  :=  est[tiae2,l]  +  gate[time2, 1] ; 

gate [time2 ,2]  :=  Abs((Sqr(mg_Bax)  -  attr[timel,obsi,l]*Sqrt(Sqr(rng_max 

-  Sqr(dmax)))/(rng_max*t2i)  -  est [time2,2] ) ; 
gate [tiae2 ,2]  :=  RMax(Abs((Sqr(rng_max)  -  attr[timei,obsi,l] 

*Sqrt (Sqr (rng_aax)  -  Sqr(dmin)))/(rng_max*t21) 

-  est [time2,2] ) ,gate[time2,2] ) ; 

gate [time2 ,2]  :=  RMax(Abs((Sqr(mg_Bin)  -  attr[timel,obsl,l] 

*Sqrt(Sqr(rng_min)  -  Sqr(dmax)))/(rng_min*t21) 

-  est[time2,2]),gate[time2,2]); 

gate [time2 ,2]  :=  RMax(Abs((Sqr(rng_min)  -  attr[timel,obsl,l] 

*Sqrt(Sqr(rng_min)  -  Sqr(dmin)))/(rng_min*t21) 

-  est[time2,2] ) ,gate[time2,2] )  +  3.0*R[2]; 
gate [time2 ,2]  :=  1 .3*gate[time2,2] ; 

{Calculate  azimuth  and  elevation  gates} 
est[time2,3]  :=  attr [timel ,obsl ,3] ; 
est [time2,4]  :  =  attr[timel,obsl,4] ; 
theta  :  =  ArcSin(dmax/est[time2,l]); 

gatertime2,3]  :=  theta/Cos(attr[timel,obsl,4])  +  3.0*R[3]; 
gate[tiBe2,4]  :=  theta  +  3.0*R[4]; 
end;  {lor  time2} 

end;  {Procedure  Calculate_SP_Estiraates} 


gate [time2 ,2] 


gate [time2 ,2]  := 


Procedure  Calculate_DP_EstiBate(timel,obsl , 

time2,obs2  :  integer; 
var  energy ,delta_energy  :  real); 

type 

vector  =  array  [1.. block]  of  real; 
var 

cluster, i, j ,k, time, atr  :  integer; 

t.obs  ••  array  [1..2]  of  integer; 

Cos_az,Cos_el,Sin_az,Sin_el,dti,r2,v2,rf  :  real; 

state, Emat,E_P  :  state.vector ; 

rmag2,rmag3,r_dot_urho  :  array  Cl.. 2]  of  real; 

mu_f actor  :  array  [1..2.1..2]  of  real; 

dt.tf  •  array  [1..3]  of  real; 

urho.rvec  •  array  Cl. .2]  of  vector; 


r_dot_part  :  array  [1..2.3..4]  o f  real; 

partial.urho  :  array  [1..2.3..4]  oi  vector; 

Jmat,J_R  :  J_«atrix; 

Pmat  :  P_matrix; 

begin 

t[l]  :=  tinel; 
t[2]  :=  time2; 
obs[l]  :=  obsl; 
obs [2]  :=  obs2; 

dt[3]  :=  frame_time[time2]  -  frame_time[timel] ; 
dti  :=  1.0/dt[3]; 
dt[2]  :=  0 . 5*dt  [3]  ; 
dt[l]  :=  -dt  [2]  ; 
tf[l]  :=  dt[l]*dti; 
tf [2]  :=  dt[2]*dti; 
for  time  :=  1  to  2  do 
begin 

Cos_az  : =  Cos (attr[t [time] , obs [time] ,3]); 

Cos_el  :=  Cos(attr[t[time] ,obs[time] ,4]) ; 

Sin_az  :=  Sin(attr[t [time] ,obs[time] ,3]); 

Sin_el  :=  Sin(attr[t[time] ,obs[time] ,4] ) ; 
urho[time, 1]  :=  Cos_el*Cos_az; 
urho[time,2]  :=  Cos_el*Sin_az; 
ur ho [time, 3]  :=  Sin_el; 
partial_urho [time ,3,1]  :=  -Cos_el*Sin_az; 
partial_urho[time,3,2]  :=  Cos_el»Cos_ar ; 
partial_urho[time,3,3]  :=  0.0; 

partial_tirho [time ,4,1]  :=  -Sin_el*Cos_az; 
p2Lrtial_urho[time,4,2]  :=  -Sin_el*Sin_az ; 
partial_urho[time,4,3]  :=  Cos_el; 
rmag2 [t ime]  : =  0.0; 
r_dot_urho [time]  :=  0.0; 
r_dot_part[time,3]  :=  0.0; 
r_dot_part [time ,4]  :=  0.0; 
for  i  :=  1  to  block  do 
begin 

rvec[time,i]  :=  attr[t[time] , obs [time] ,l]*ur ho [time, i]  +  Sxs [t [time] , i] ; 
rmag2[time]  :=  rmag2[time]  +  Sqr (rvec [time , i] ) ; 

r_dot_urho[time]  :=  r_dot_urho[time]  +  rvec [time, i]*urho[t ime, i] ; 
for  atr  :=  3  to  4  do 

r_dot_part[time,atr]  :=  r_dot_part [time, atr] 

+  rvec[time,i]*partial_urho[time,atr, i] ; 

end;  {for  i> 

rmag3[time]  :=  Sqrt (rmag2 [time] )*rmag2 [time] ; 
mu_f actor [time, 1]  :=  mu*dt[l]*dt[2]/(2.0*rmag3[time]) ; 
mu_factor[time,2]  :=  3. 0*mu _f actor [t ime, l]/rmag2 [time] ; 
mu_f actor [time, 1]  :=  1.0  -  mu_f actor [time, 1] ; 
end;  {for  time} 
r2  :=  0.0;  v2  :=  0.0; 


for  i  :=  1  to  block  do 
begin 

j  :=  i  +  block; 

state[i]  :=  tf [2]*mu_factor[i,l]*rvec[l,i] 

-  tf  [l]*mu_factor[2,l]*rvec[2,i] ; 
r2  :=  r2  +  Sqr(atate[i] ) ; 
state [j]  :=  dti*(mu_factor[2,l]*rvec[2,i] 

-  mu_f actor [l,l]*rvec[l,i]); 
v2  :=  v2  +  Sqr(atat«[j]) ; 

Jmat[i,l]  :=  tf [2]*(mu_factor[l,l]*urho[l,i] 

+  mu_f actor [1 , 2] *r_dot_urho  Cl] *rvec [1 , i] ) ; 

Jmat [i, 2]  :=  0.0; 

Jmat[i,6]  :=  -tf [l]*(mu_factor[2,l]*urho[2,i] 

+  mu_factor[2,2]*r_dot_urho[2]*rvec[2,i]) ; 

Jmat  [i ,  6]  :  =  0 . 0 ; 

Jmat[j,l]  :=  -dti*(mu_factor[l,l]*urho[l,i] 

+  mu_factor[l,2]*r_dot_urho[l]*rvec[i,i]) ; 

Jmat  [  j  ,  2]  :  =  0 . 0 ; 

Jmat [j, 6]  :=  dti*(mu_factor[2,l]*nrho[2,i] 

+  mu_f  actor  [2 , 2]  *r_dot_urho  [2]  *rvec  [2 ,  i]  ) ; 

Jmat  C j , 6]  : =  0 . 0 ; 
for  atr  :=  3  to  4  do 
begin 

Jmat [i, atr]  :=  attr Ctimel ,obsl ,l]*tf [2] 

*  (mu_f  actor  [1,1]  *part ial_urho  [1 ,  atr ,  i] 

+  mu.factorCl^Jer.dot.partCl.atri^rvecCl.i]) ; 
Jmat[i,atr+4]  :=  -attr [time2,obs2, l]*tf [1] 

* (mu_f actor [2,1] +partial_urho [2 , atr , i] 

+  mu_f actor [2,2] *r_dot_part [2 ,atr] +rvec [2 , i] ) ; 
Jmat[j,atr]  :=  -attr[timel,obsl,l]*dti 

*(mu_f actor [1 , 1] *partial_urho [1 , atr , i] 

+  mu_f actor [1,2] *r_dot_part [1 , atr] *rvec [1 , i] ) ; 
Jmat[j ,atr+4]  :=  attr[time2,obs2,l]*dti 

* (mu_f actor [2 , 1] *part ial.urho [2 , atr , i] 

+  mu_factor[2,2]*r_dot_part[2,atr]*rvec[2,i]) ; 

end;  {for  atr} 
end;  {for  i} 

energy  :=  v2/2.0  -  mu/Sqrt(r2); 
rf  :=  mu/(r2*Sqrt(r2)) ; 
for  i  :=  1  to  block  do 
begin 

j  :=  i  +  block; 

EmatCi]  :=  state [i]*rf; 

Emat[j]  :=  state[j]; 
end;  {for  i} 
for  i  :=  1  to  nest  do 
for  j  :=  1  to  nterms  do 
begin 
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lor  k  :=  1  to  ntenns  do 

:=  +  Jmat[i,k]*Rvar[2,k, j] ; 

end;  {lor  j> 
lor  i  :=  1  to  nest  do 
lor  j  :=  1  to  nest  do 
begin 

PnatCi.j]  :=  0.0; 
lor  k  :=  1  to  ntenns  do 

PmatCi.j]  :=  P»at[i,j]  +  J_R[i,k]*Jmat[j  ,k]  ; 
end;  {lor  j> 
lor  j  :=  1  to  nest  do 
begin 

EJP[j3  :=  0.0; 

lor  k  :=  1  to  nest  do 

E_P  Cj]  :=  E_P  [j]  +  Eaat  [k]  *Pmat  [k ,  j]  ; 
end;  {lor  j> 
delta_energy  :=  0.0; 
lor  j  :=  1  to  nest  do 

delta_energy  :=  delta_energy  +  E_P[j]*Emat[j] ; 
end;  {Procedure  Calculate_DP_Estimate} 


cluster, i, j ,k,ti»e,atr 
t  ,obs 

Cos.az , Cos_el , Sin_az , Sin_el ,r2 , v2 ,rl 
Emat,E_P 

mag  2 ,  mag  3 ,  r  _dot  _urho 

r_dot_part 

mu_lactor 

drho 

dt.tl 

urho.rvec 

partial.urho 

Jmat, J_R 


:  integer; 

:  array  [1..3]  ol  integer; 

:  real; 

:  state_vector ; 

:  array  [1..3]  ol  real; 

:  array  [1..3.3..4]  ol  real; 

:  array  [1..3.1..2]  ol  real; 

:  vector; 

:  array  [1..3]  ol  real; 

:  array  [1..3]  ol  vector; 

:  array  [1..3.3..4]  ol  vector; 
:  J_matrix; 


begin 

t[l]  :=  timel;  t[2]  :=  time2;  t[3]  :=  time3; 
obs[l]  :=  obsl;  obs[2]  :=  obs2;  obs[3]  :=  obs3; 
dt[l]  :=  lrane_time[time2]  -  lrajne_tiBe[timel] ; 
dt[2]  :=  lrame_time[time3]  -  lraae_time [time2] ; 
dt[3]  :=  dt  [1]  +  dt  [2]  ; 
tl[l]  :=  -dt[2]/(dt[l]*dt[3]) ; 


I  j.l'lJi’i.H-l  4.*  t.»  »,*  l,«  h* JiLk1  I 


tl  [2]  :=  (dt  [2]-dt  [1]  )/ (dt  [1]  *dt  [2])  ; 
tl  [3]  :=  dt  [1]  /  (dt  [2]  ♦dt  [3]  ) ; 
lor  time  : =  1  to  3  do 
begin 

Cos.az  Cos (attrCt [time] ,obs[time] ,3] ) ; 

Cos_el  :=  Cos (attr[t [time] ,obs[time] ,4]) ; 

Sin.az  :=  Sin(attr[t[time] ,obs[time] ,3]); 

Sin_el  :=  Sin(attr[t [time] ,obs [time] ,4] ) ; 
urho[time,l]  :=  Cos_el*Cos_az; 
urho[time,2]  :=  Cos_eleSin_az; 
urho[time,3]  :=  Sin_el; 
partial_urho [time, 3,1]  :=  -Cos_el*Sin_az; 
partial_urho[time,3,2]  :=  Cos_el*Cos_az; 
partial_urho[time,3.3]  :=  0.0; 

partial_nrho [time, 4,1]  :=  -Sin_el*Cos_az; 
partial_urho[time,4,2]  :=  -Sin_el*Sin_az ; 
partial_urho[time,4,3]  :=  Cos_el; 
end;  {lor  time] 
r2  :=  0.0; 
v2  :=  0.0; 

lor  i  :=  1  to  block  do 
begin 

j  :=  i  +  block; 

state [i]  :=  attr[time2,oba2,l]*urho[2,i]  +  Sxs [time2, i] ; 
r2  :=  r2  +  Sqr(state[i] ) ; 

drho[i]  :=  tl[l]*urho[l,i]  +  tl [2] ♦urho [2 , i]  +  tl[3]*urho[3,i] ; 
state [j]  :*  attr[time2,obs2,2]*urho[2,i]  +  attr[time2,obs2,l]*drho[i] 
+  Sxs[time2, j] ; 
v2  :=  v2  +  Sqr(8tate[j] ) ; 

Jmat  [i ,  1]  :  *  urho  [2 ,  i]  ; 

Jmat[i,2]  :=  0.0; 

Jmat [i , 3]  : =  0 . 0 ; 

Jmat[i,4]  :=  attr[time2,obs2,l]*partial_urho[2,3,i] ; 

Jmat [i, 5]  :=  0.0; 

Jmat  [i ,  6]  :  =  0 . 0 ; 

Jmat[i,7]  :=  attr[time2,obs2,l]*partial_urho[2,4,i] ; 


Jmat [i , 8] 
Jmat  [j ,  1] 
Jmat  [j  ,2] 


0.0; 

drho[i] ; 
urho  [2 ,  i]  ; 


lor  atr  :=  3  to  4  do 
begin 

k  :=  3* (atr  -  2) ; 

Jmat [j ,k]  :=  attr[time2,obs2,l]*tl[l]*partial_urho[l,atr,i] ; 

Jmat[j,k+1]  :=  (attr[time2,obs2,2]  +  attr[time2,obs2,l]*tl[2]) 
♦partial_urho[2,atr,i]  ; 

Jmat[j,k+2]  :=  attr[time2,obs2,l]*tl [3]*partial_urho[3,atr,i] ; 


Jmat[j,k+2]  := 
end;  {lor  atr} 
end;  {lor  i} 
energy  :=  v2/2.0 


mu/Sqrt(r2) ; 
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rf  :=  mu/(r2*Sqrt(r2)) ; 
for  i  :=  1  to  block  do 
begin 

j  :=  i  +  block; 

EnatCi]  :=  state [i]*rf; 

EmatCj]  :=  state Cj]; 
end;  {for  i} 

{Fora  state  covariance  matrix} 
for  i  :=  1  to  nest  do 
for  j  :=  1  to  nterms  do 
begin 

•J_R[i,j]  :=  0.0; 

for  k  :=  1  to  nterms  do 

J_R[i,j]  J_R[i.j]  +  Jmat[i,k]*Rvar [3,k, j] ; 
end;  {for  j} 
for  i  :=  1  to  nest  do 
for  j  :=  1  to  nest  do 
begin 

Pmat[i,j]  :=  0.0; 
for  k  1  to  nterms  do 

Pmat[i,j]  :=  Pmat[i,j]  +  J_R[i,k]*Jmat[j ,k] ; 
end;  {for  j} 
for  j  :=  1  to  nest  do 
begin 

E_P[j]  :=  0.0; 

for  k  :=  1  to  nest  do 

E_P[j]  :=  E_P[j]  +  Emat [k] *Pmat [k ,  j]  ; 
end;  {for  j> 
delta.energy  :=  0.0; 
for  j  :=  1  to  nest  do 

delta_energy  :=  delta_energy  +  E_P[j]*Emat[j] ; 
end;  {Procedure  Calculate_TP_Estimate> 


Procedure  Perform_Cluster_Initiation; 
label  1; 
type 

pairs  =  record 

timel,obsl,time2,obs2  :  integer; 
end;  {record} 
triples  -  record 

timel>obsl,tiae2,obs2,time3,obs3  :  integer; 


metric 


:  real; 


end;  {record} 


boolean; 


time,obs,atr,ntime, start .stop, 
count, ocountl ,ocount2, cluster, 
missed.gates , po inti, point2, arc  :  integer; 
delta, specif ic.energy  :  real; 


active 

estimates, gates 
index , back.index , incident 
pointer 
pair 
triple 

Procedure  Solve.QP (sores 
var  best 

integer; 
real; 
bvector; 

boolean; 


:  bvector; 

:  0 .matrix; 

:  array  [spsui]  of  ivector; 

:  array  [framel.. -1]  of  ivector; 

:  array  [framel. .-1,1. .clusters]  of  pairs; 
:  array  [1 . .clusters]  of  triples; 

integer; 
bvector) ; 

var 
i 

■incost 
solution 

Function  Feasible 
label  1; 
var 

result  :  boolean; 
time, obs  :  integer; 
begin 

result  :=  true; 
for  time  :=  framel  to  0  do 

for  obs  :=  1  to  nr_unassigned[time]  do 
if  incident [time, obs]  >  1  then 
begin 

result  : =  false; 
goto  1; 
end;  {if> 

Feasible  :=  result; 
end;  {Function  Feasible] 

Procedure  Search (j  :  integer; 

c  :  real) ; 

var 

k  :  integer; 
d  :  real ; 
begin 

for  k  :=  j+1  to  arcs  do 
begin 

solution Dt]  :=  false; 
with  triple [k]  do 
begin 

Decrement (incident [time 1 ,back_index[timel ,obsl]] ) ; 

Decrement (incident [time2 , back. index [time2 , obs2] ] ) ; 

Decrement (incident [time3 , back.index [time3 , obs3] ] ) ; 
d  :=  c  -  metric; 
end;  {with} 
if  not  Feasible  then 
SeaLrch(k,d) 
else 

if  d  <  mincost  then 
begin 
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best  :=  solution; 
nincost  :=  d; 
end;  {if} 

solution Dc]  :  =  true; 
with  triple [k]  do 
begin 

Increment ( incident [t ime 1 , back. index [t imel , obs 1] ] ) 
Increment ( incident [time2 , back.index [t ime2 , obs2] ] ) 
Increment (incident [t ime 3 , back.index [t ime3 , obs 3] ] ) ; 
end;  {with} 
end;  {for} 

end;  {Procedure  Search} 
begin 

solution  :=  best; 
mincost  :=  big; 
if  not  Feasible  then 
Search(O.O) ; 

end;  {Procedure  Solve.QP} 
begin 

Start_Timer(5) ; 

for  time  :=  framel  to  0  do  {Index  unassigned  observations} 
begin 
obs  : =  0 ; 

for  count  :=  1  to  nr .unassigned [time]  do 
begin 
repeat 

Increment (obs) ; 
back.index [time, obs]  :=  0; 
until  not  as signed [time, obs ] ; 
index [time , count]  : =  obs ; 
back.index [time, obs]  :=  count; 
end;  {for  count} 
end;  {for  time} 

for  staxt  :=  framel  to  -1  do  {Determine  possible  pairr" 
begin 
obs  : =  0 ; 

for  ocountl  :=  1  to  nr .unassigned [start]  do 
begin 

pointer[start .ocountl]  :=  obs  +  1; 
Calculate_SP_Estimates( start .index [start .ocountl] , 

estimates, gates) ; 
for  stop  :=  start+1  to  0  do 
begin 

for  ocount2  :=  1  to  nr .unassigned [stop]  do 
begin 

missed_gates  :=  0; 
for  atr  :=  1  to  attributes  do 
begin 


delta  :  =  Abs (estimates [atop, atr] 

-  attr [stop , index [stop , ocount2] , atr] ) ; 
if  delta  >  gates [stop, atr]  then 
if  (delta  >  2.0*gates[stop,atr] ) 
or  (missed_gates  >=  max.missed)  then 
goto  1 
else 

Increment (miss ed_gates) ; 
end;  {for  atr} 

Calculate_DP_Estimate (start , index [start , ocountl] , 

atop .index [stop ,ocount2] , 
specif ic_energy .delta) ; 

if  specif ic_energy  <=  max_energy  +  3.0*Sqrt(delta)  then 
begin 

Increment (obs) ; 
frith  pair  [start,  obs]  do 
begin 


timel  :=  start; 

obsl  :=  index [start .ocountl] ; 
time2  :=  stop; 

obs2  :=  index [stop, ocount 2] ; 
end;  {with} 
end;  {if} 

1:  end;  {for  ocount2} 

end;  {for  stop} 
end;  {for  ocountl} 

pointer  [start, nr ..unassigned [start] +1]  :=  obs  +  1; 
pointer  [start ,0]  :=  obs; 
end;  {for  start} 

{Form  all  possible  triples  and  calculate  metric;  eliminate  infeasible  triples} 
for  time  :=  frame 1  to  0  do 

for  obs  :=  1  to  nr_unassigned[time]  do 
incident [time, obs]  :=  0; 
count  : =  0 ; 

cluster  :=  nr_clusters [0] ; 
for  start  :=  framel  to  -2  do 
begin 

for  ocountl  :=  1  to  pointer [start.O]  do 
begin 

stop  : =  pair [start, ocount 1] . time2; 
if  stop  <>  0  then 
begin 

obs  :=  back_index [stop, pair [start .ocountl] .obs2] ; 
pointl  :=  pointer [stop, obs] ; 
point2  :=  pointer[stop,obs+l] ; 
for  ocount2  :=  pointl  to  point2-l  do 
begin 

Increment (cluster) ; 
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Calculat « _TP_Est imat e ( start , pair [start , ocount 1] . obs 1 , 

stop, pair  [stop, ocount 2] . obsl, 
pair[stop,ocount2] .time2, 
pair [stop, ocount 2] .obs2, 

Sx [cluster] ,P [cluster] , 
specif ic.energy , delta) ; 

if  specif ic_energy  <=  max_energy  +  3.0*Sqrt(delta)  then 
begin 

Increment (count) ; 
nr .missing [cluster]  :=  0; 
with  triple [count]  do 
begin 

timel  :=  start; 

obsl  :=  pair [start, ocount l] .obsl; 
time2  : =  stop; 

obs2  :=  pair [stop , ocount 2] . obsl ; 
time3  :=  pair [stop, ocount 2] .time2; 
obs3  :=  pair [stop, ocount2] .obs2; 
metric  :=  specif ic.energy ; 

Increment (incident [timel ,back_index [timel .obsl]] ) ; 
Increment(incident[time2,back_index[time2,obs2]] ) ; 
Increment ( incident [t ime3 , back.index [time3 , obs3] ] ) ; 
active [count]  :=  true; 
end;  {with} 
end;  {if} 

end;  {for  ocount2> 
end;  {if} 

end;  {for  ocountl} 
end;  {for  start} 

{Solve  remaining  quadratic  program  using  implicit  enumeration} 
Solve_QP(count, active) ; 
cluster  :=  nr .clusters [0] ; 
move  :=  false; 
for  arc  :=  1  to  count  do 
if  active[arc]  then 
begin 

Increment (cluster) ; 
nr .missing [cluster]  :=  0; 
if  move  then 
begin 

Sx[cluster]  :=  Sx[nr_clusters[~l]+arc] ; 

P[cluster]  :=  P [nr .cluster s [- 1] +arc] ; 
end;  {if  move} 
with  triple[arc]  do 
begin 

as signed [timel ,obsl]  :=  true; 
Decrement(nr_unassigned[timel] ) ; 
observation  [timel  .duster]  :=  obsl; 
assigned [time2 ,obs2]  :=  true; 
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Decrament(nr_unas8igied[time2] ) ; 
observation [time2, cluster]  :=  obs2; 
assigned[time3,obs3]  :  =  true; 

Decrement (nr .unass igned[time3] ) ; 
observation[time3, cluster]  :=  obs3; 
for  time  :  =  timel  to  0  do 
begin 

status [time .cluster]  :=  time2  -  1; 

Increment (nr .clusters [time]) ; 

Increment (nr .active [time] ) ; 

if  (time  <>  timel)  and  (time  <>  time2)  and  (time  <>  time3)  then 
observation [time, cluster]  :=  0; 
end;  {for  time) 
end;  {with} 
end  {if} 
else 

move  : =  true ; 

Stop.Timer(S) ;  {Ti: 

end;  {Procedure  Perf orm.Cluster .Initiation} 


{Timing} 


{***  Outputs  ****************************************************************} 


Procedure  Echo_Cluster_Assignment(time  :  integer); 


cluster .count  :  integer; 
begin 

if  frame_time[time]  >=0.0  then 
begin 

srite(tcldata,frame_time[time] :7: 1) ; 
cluster  :=  0; 

for  count  :=  1  to  nr .active [time]  do 
repeat 

Increment(cluster) ; 

if  status [time, cluster]  <  0  then 

Vr it  e ( t cldata , target  s [t ime , oba  ervat ion [t ime , clust  er] ] : 4) 
else 

Write(tcldata, ’  ’); 

until  status [time, cluster]  <  0; 

Writeln(tcldata) ; 
end;  {if} 

end;  {Procedure  Echo.Cluster.Assignment} 


Procedure  Output .Cluster.Residuals ; 


count, cluster, target, i,j  :  integer; 

position, velocity  :  real; 

Sxt  :  array  [1 .. clusters]  of  state.vector ; 
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begin 

{Input  target  states  for  current  observation  frame} 
repeat 

Sxt [current .target. number]  :=  current .target. values; 

Read(tsdata, current .target) ; 
until  (current .target. time  >  frame_time[0] ) ; 

{Compare  estimated  to  true  target  states} 
cluster  : =  0; 

for  count  :=  1  to  nr .active [0]  do 
repeat 

Increment(cluster) ; 
if  status [time .cluster]  <  0  then 
begin 

target  :=  targets[0, observation^, cluster]]  ; 
position  :=  0.0; 
velocity  :=  0.0; 
if  target  >  0  then 
begin 

for  i  :=  1  to  block  do 
begin 

j  :=  i  ♦  block; 

position  :=  position  +  Sqr(Sx[cluster , i]  -  Sxt [target , i] ) ; 
velocity  :=  velocity  +  Sqr(Sx[cluster , j]  -  Sxt[target, j]) ; 
end;  {for  i} 

position  :=  Sqrt (position) ; 
velocity  :=  Sqrt (velocity) ; 
end;  {if  target  >  0} 

Vrite(c8data, cluster : 4, target :4, frame .time [0] :7: 1) ; 

Wr iteln(csdata, posit ion: 11: 1, velocity :7: 1) ; 
end;  {if  status  <  0} 
until  status [0, cluster]  <  0; 
end;  {Procedure  Output .Cluster .Residuals} 

{***  Main  Program  **********+****+************************+******************} 

BEGIN 

{*****************************} 

{**  Perform  Initializations  **} 

{*****************************} 

Init .Times; 

Start _Timer(0) ; 

Init .Program; 

Initialize_RK78; 

Init ialize.Estimat ion; 

Initialize.Clustering; 


{Timing} 

{Timing} 


{***********+***********************} 

{**  Perlorm  sequential  clustering  **} 

{***********************************} 

repeat 

Input _Data; 

Forecast; 

Calculate_Metrics ; 

Perlorm.Cluster .Assignment; 

Update .Estimates ; 

Perform.Cluster .Initiation; 

Echo.Cluster.AssignmentCframel) ; 

Output .Cluster .Residuals ; 
until  EOI; 

■(*♦*****************************> 

{**  End  sequential  clustering  **} 

{*******************************} 
lor  time  :=  lramel+1  to  0  do 
Echo .Cluster .Assignment(time) ; 

Stop.Timer(O) ;  {Timing} 

Report_Times(6) ;  {Timing} 
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